Surface charge relaxation controls the lifetime of out-of-equilibrium colloidal crystals

Abstract

Interactions between charged colloidal particles are profoundly influenced by charge regulation and charge renormalization, rendering the effective potential highly sensitive to local particle density. In this work, we investigate how a dynamically evolving, density-dependent Yukawa interaction affects the stability of out-of-equilibrium colloidal structures. Motivated by a series of experiments where unexpectedly long-lived colloidal crystals have suggested the presence of like-charged attractions, we systematically explore the role of charge regulation and charge renormalization. Using Poisson–Boltzmann cell theory, we compute the effective colloidal charge and screening length as a function of packing fraction. These results are subsequently incorporated into Brownian dynamics simulations that dynamically resolve the evolving colloid charge as a function of time and local density. In the case of slow relaxation dynamics, our results show that incorporating these charging effects significantly prolongs the lifetimes of out-of-equilibrium colloidal crystals, providing an explanation for the experimental observation of long-lived crystals. These findings demonstrate that the interplay of surface charge dynamics and colloidal interactions can give rise to complex and rich nonequilibrium behavior in charged colloidal suspensions, opening new pathways for tuning colloidal stability through electrostatic feedback mechanisms.

---

1 Introduction

Suspensions of charge-stabilized colloids have attracted sustained interest in soft matter research for nearly a century, driven by their relevance in diverse applications ranging from industrial coatings to pharmaceutical formulations.¹ The theoretical framework for understanding interactions in these systems is provided by the Derjaguin–Landau–Verwey–Overbeek (DLVO) theory, which describes the screened electrostatic colloidal interactions in an electrolyte through a distance-dependent Yukawa potential.^2,3 Despite its conceptual simplicity, the DLVO potential has proven highly effective in describing the behavior of most charged colloidal systems and remains a central tool in colloid science.^4,5

Since the 1980s, several experiments have cast doubts on the validity of DLVO theory for colloidal suspensions under low-salt conditions.^6–8 These studies reported evidence of long-ranged like-charge attractions between colloids—an effect not predicted by the DLVO framework. Observed phenomena suggest such attractions include vapor–liquid phase separation,⁹ the formation of dilute voids,^10–12 clustering of colloids,^13,14 and colloidal crystals with unexpectedly long lifetimes.^15,16 The origin of these apparent attractive interactions has yet to be resolved. It is important to note, however, that these findings remain controversial. Some results have proven difficult to reproduce,¹⁷ while others may stem from experimental artifacts, such as out-of-equilibrium hydrodynamic effects caused by the proximity of the particles to a substrate.¹⁸

Another intriguing phenomenon observed in charged-colloid suspensions is reentrant melting, where colloids undergo a phase transition from a crystal to a fluid phase upon increasing the density.^19,20 This behavior seems to originate from charge regulation mechanisms. Here, charge regulation is an umbrella term that includes a broad class of processes in which neither the surface charge nor the surface potential is constant.^21,22 The specific form of charge regulation depends on the underlying charging mechanism, which can involve either ionizable surface groups²² or the adsorption of ionic species.^19,20,23 Depending on the nature of the surface chemistry and solvent environment, the surface charge may increase²⁴ or decrease^19,20,23,25 with increasing colloid density. Furthermore, these charging mechanisms do not occur instantaneously, but are typically dynamic with finite timescales associated with charging processes.^26–29

A related but distinct density-dependent phenomenon is charge renormalization, in which the bare colloid charge is replaced by a reduced effective charge that accounts for the counterions condensed on the particle surface.^30–32 This effective charge enables the use of linearized Poisson–Boltzmann theories, while still capturing nonlinear screening effects. As a result, both the colloid charge and the screening length become renormalized quantities.^30,31 Charge renormalization is widely used to reconcile the significantly lower effective charges observed in experiments with the much higher bare charges typically obtained from titration measurements.^33–37

In this work, we employ Poisson–Boltzmann cell model calculations to determine the renormalized charge and screening length as a function of the packing fraction of the colloids,^23,25,31 following the approach of Zoetekouw et al.²⁵ These results serve as input for Brownian dynamics (BD) simulations that explicitly account for the time- and density-dependent evolution of the colloid charge. Using these simulations, we investigate the structural lifetimes of colloidal crystals in both cubic and planar slab geometries, focusing on how they are influenced by variations in charging dynamics.

To study this theoretical model in a physically relevant setting, we chose parameters consistent with the experiments of Larsen and Grier.^15,16 In their experiments, polystyrene sulfate colloids with a diameter of 652 ± 5 nm were used at a low volume fraction of η = 0.02.¹⁶ These particles carried a high surface charge due to strongly acidic sulfate surface groups, enabling them to be compressed into crystalline structures using an oscillating electric field. Upon removal of the field, the crystals reverted to a homogeneous fluid phase. Notably, superheated colloidal crystals prepared at relatively high ionic strength melted within approximately 10 seconds, whereas those at lower ionic strength remained metastable for up to 30 minutes, before reaching their equilibrium state.^15,16 These long-lived structures also exhibited significant density variations, with differences as high as 70%.¹⁵

We investigate these so-called superheated crystals from a new perspective, aiming to demonstrate not only our new method for incorporating dynamic surface charge relaxation but also the impact of density-dependent interactions on the stability of charged colloidal crystals in suspension. Using the BD simulations described above, we systematically study how the lifetime of a crystal depends on charge regulation and charge renormalization, with particular focus on the rate at which these charging and decharging processes occur.

We demonstrate that crystal lifetimes are significantly enhanced by reducing the rate at which colloids equilibrate their surface charge, indicating that slow, density-dependent charging processes contribute to stabilizing the superheated crystal state. However, the simulated crystals do not exhibit the pronounced density differences as observed in experiments,¹⁵ suggesting that additional mechanisms beyond our current model are needed to fully explain the experimental observations.

This paper is organized as follows. In Section 2, we introduce the experimental and theoretical background motivating our Poisson–Boltzmann (PB) cell theory calculations and discuss the physically relevant parameter regimes. Section 3 outlines the BD simulations and the machine-learning-based analysis techniques employed in this work. In Section 4, we present the results of our PB cell theory calculations, which are then employed in BD simulations to investigate the structural lifetimes of superheated colloidal crystals in both cubic and planar slab geometries, with the latter compared directly to experimental observations. Finally, we summarize our key findings and conclusions in Section 5.

2 Model and theory

2.1 Interactions between charged colloids

According to the classical DLVO theory, the effective interaction potential between two identical, homogeneously charged colloidal spheres of radius a and total charge Ze with e the elementary charge, immersed in a solvent containing both coions and counterions with dielectric constant ε and Debye screening length κ⁻¹, is described by a purely repulsive screened-Coulomb (Yukawa) interaction^2,3where r represents the center-of-mass distance between the spheres, β = 1/k_BT with k_B Boltzmann's constant and T the temperature, and λ_B = e²/4πεk_BT is the Bjerrum length. Note that we ignore the dispersion forces here.

In addition, we describe the hard-core repulsion between the colloids by a Weeks–Chandler–Andersen (WCA) potential³⁸

where σ = 2a denotes the particle diameter with a the particle radius, and ε_WCA is a parameter that sets the strength of the WCA potential. We use βε_WCA = 40, which has been used extensively in previous studies to model hard spheres.³⁹ The total interaction potential between two colloids is simply the sum of the two potentials, βU(r) = βU_Yuk + βU_WCA(r). In this work, the range of the colloid interaction always exceeds r > 2^1/6σ, so the specific form of the WCA potential serves primarily as a formal repulsive core rather than playing a significant physical role.

This study explores the effect of charge renormalization and charge regulation on the stability of out-of-equilibrium colloidal structures by taking into account a dynamically evolving Yukawa interaction that depends on the instantaneous local colloid density. Below, we describe how charge renormalization and charge regulation are accounted for within the framework of a Poisson–Boltzmann cell theory.

2.2 Charge renormalization and charge regulation

Originally, Derjaguin, Landau, Verwey, and Overbeek derived a Yukawa-like interaction potential by linearizing the Poisson–Boltzmann equation and applying the linear superposition approximation, i.e., using single-sphere solutions of the linearized Poisson–Boltzmann (PB) equation to describe the interaction potential between two spheres. This provided a convenient analytical form for screened electrostatic interactions between charged colloids.^2,3 However, when the electric or surface potential of the colloids become sufficiently large, non-linear effects become significant and the linearization of the Poisson–Boltzmann equation and the linear superposition break down. In such regimes, experimental studies have consistently shown that the effective colloid charge, obtained by fitting the interactions to a Yukawa potential, is significantly lower than the bare charge determined through titration, which directly probes the surface chemical groups.^33–35 In other words, the bare charge overestimates the charge needed in the Yukawa potential to accurately describe phase behavior. To resolve this, Alexander et al. introduced the concept of charge renormalization, in which the long-range asymptotic decay of the full non-linear PB solution is matched to a Yukawa form with an effective charge.^30,31 This approach captures non-linear screening effects while retaining the analytical simplicity of the Yukawa potential, thereby improving the accuracy of theoretical predictions for colloidal interactions and phase behavior.

In many theoretical treatments, the charge on a colloidal surface is assumed to be constant and fixed. However, this assumption often fails to capture the behavior observed in realistic experimental systems, where the origin of surface charge must be explicitly considered. Typically, the charge arises from the dissociation of chemical groups on the colloidal surface in a polar solvent, such as water, resulting in a net surface charge. As a result, the colloid charge is not a fixed quantity but a dynamic one that can vary with system parameters such as temperature, pH, colloid packing fraction, and background ion concentration. This phenomenon, known as charge regulation, must be taken into account to accurately describe colloidal interactions under realistic experimental conditions.

2.2.1 Poisson–Boltzmann cell theory. To account for non-linear screening and charge regulation effects, we employ the spherical cell approximation, which allows us to solve the full non-linear Poisson–Boltzmann (PB) equation. In this framework, a system of N colloidal particles is modeled by partitioning the volume into identically sized spherical Wigner–Seitz cells, each containing a single colloid at its center.⁴⁰ The cell is characterized by a radius R, which is related to the particle size a and the colloid packing fraction η through the relation η = η_max(a/R)³. Here, η_max ≃ 0.74 represents the maximum packing fraction of a face-centered cubic (FCC) or hexagonal close-packed crystal of hard spheres, as used in ref. 25. While this convention of introducing a maximum packing fraction is not universally adopted in the literature,^23,31,41 it can be reverted to the standard form by setting η_max = 1. This reduces the complex many-body problem to a one-body problem with spherical symmetry, enabling the numerical solution of the non-linearized PB equation inside a single spherical cell. Within the Poisson–Boltzmann framework, the ion distributions follow a Boltzmann distribution n_±(r) = n^∞ exp(∓ϕ(r)), where ϕ(r) = eψ(r)/k_BT is the dimensionless electric potential, r is the radial distance from the center of the colloid, and n^∞ is the bulk concentration of coions or counterions.⁴² Substituting the total charge density ρ_c(r) = en₊(r) − en₋(r) into the Poisson equation from electrostatics ∇²ψ(r) = −ρ_c(r)/ε, leads to the non-linear PB equation

∇²ϕ(r) = κ² sinh ϕ(r), (3)

where is the inverse Debye screening length.^42,43

To determine the renormalized charge, we follow the approach by Alexander et al.³⁰ and numerically solve the non-linear Poisson–Boltzmann (PB) eqn (3) under appropriate boundary conditions. Since the PB equation is a second-order differential equation, two boundary conditions are required to obtain a unique solution. The first boundary condition,

ensures total charge neutrality within the spherical cell.⁴³ The second boundary condition is imposed at the colloidal surface and reflects the surface chemistry of the colloids. This condition can vary depending on whether the surface charge or surface potential is held fixed, or if charge regulation is taken into account. In this work, we focus on dissociation reactions of acidic groups, in line with the polystyrene sulfate colloids we use as a model system.^15,16 This dissociation equilibrium is described by

AH ⇌ A⁻ + H⁺, (5)

where A⁻ represents a surface charge group on the colloid and H⁺ the released counterion. This reaction leads to a relation between surface charge and surface potential^22,44where M is the number of surface groups, Z is the bare surface charge and pK_a and pH = −log([H⁺]) are the negative logarithmic acid dissociation constant and the negative logarithmic hydrogen ion concentration, respectively. The derivation of this relation, presented in Appendix A, follows directly from the chemical equilibrium conditions at the colloid surface. The bare charge is related to the electric potential via Gauss’ law,^41,43 forming the second boundary condition to the non-linear PB equation. For strongly acidic surface groups, the assumption of a constant charge is equally viable. Moreover, ion-exchange resins are used to replace stray ions with a dominant species. Thus, for all intents and purposes, the counterion concentration can be taken to be equal to [H⁺].

2.2.2 Linearized Poisson–Boltzmann cell theory. Using the solution to the non-linear Poisson–Boltzmann equation, we can calculate the renormalized charge Z* and the renormalized inverse screening length . However, when transitioning from the non-linear to the linearized Poisson–Boltzmann equation, one must choose a reference point about which the equation is linearized. Conventionally, this is done at the cell boundary R,^23,41 but one could also linearize around the mean electric potential ϕ_m, resulting in ref. 31 and 45andwhere the mean electric potential ϕ_m is given byHere, V_f refers to the free volume within the spherical Wigner–Seitz cell, which reads^31,45In this work, we employ the cell model linearized around the mean electric potential (eqn (7)), as it enables a one-to-one mapping between Z* and over a broad range of packing fractions η. This mapping simplifies the numerical methods concerned with updating the equations of motion for these two parameters, as will be explained in more detail in Section 3.1. The theoretical framework and results for the linearized cell model using the cell boundary R are provided in Appendix B.

2.3 Relaxation dynamics of surface charge and screening length

Using eqn (7) and (8), the renormalized surface charge Z* and inverse screening length are evaluated as functions of the local colloid packing fraction η. These quantities form the basis for the out-of-equilibrium Brownian dynamics simulations, which employ the Yukawa potential defined in eqn (1). To determine the local colloid packing fraction η, the simulation volume is partitioned into Wigner–Seitz cells using the Voronoi tesselation, for instance via voro++.⁴⁶ As the system evolves, the local packing fraction η—and consequently the surface charge Z* and reduced screening parameter —changes dynamically. Instead of updating these quantities instantaneously, we introduce a time delay that accounts for the finite relaxation time between changes in the chemical environment and the relaxation to a new (renormalized and regulated) charge. This rate of change is governed by the surface chemistry underlying charge renormalization and charge regulation.

2.3.1 A microscopic model for charge renormalization. So far, we have treated charge renormalization in a coarse-grained manner. To develop a more microscopic description, we draw parallels to Manning's counterion condensation theory for polyions.⁴⁷

Manning's 1969 theory predicts that counterions condense onto an infinitely long, thin line charge when the spacing between charged groups is smaller than the Bjerrum length.⁴⁷ A later refinement classified condensed counterions as free, territorially bound, or site-bound.⁴⁸ Territorially bound ions can move freely along the polyion surface, whereas site-bound ions are more restricted. However, experiments and simulations show that the site binding is short-lived, typically less than a few nanoseconds.^49–51 For this reason, we only consider territorially bound and free counterions.

Manning condensation was formulated for cylindrical polyions, and counterion behavior depends strongly on geometry.^32,52–54 This geometric dependence is well illustrated by toy models such as those of Tang and Rubinstein.⁵³ Despite this dependence on geometry, counterion condensation and charge renormalization are closely related.^33,54–56 For example, Diehl and Levin⁵⁵ used molecular dynamics simulations with monovalent counterions around a spherical colloid and recovered the renormalized charge as predicted by Alexander et al.³⁰ In their approach, condensation was identified by comparing the potential energies and kinetic energies of the counterions, and the renormalized charge was obtained by substracting the number of condensed counterions from the bare charge.

2.3.2 Adsorption kinetics for counterion condensation. In light of these observations, we propose that charge renormalization can be interpreted as an adsorption process. This interpretation is consistent with its established role in models of dynamical charge regulation.^57–59 Adsorption-based frameworks are also widely used in studies of linear flexible polyions,^50,60–63 and, to a lesser extent for globular polyions.⁶⁴

A wide range of models exists for describing adsorption kinetics, varying in complexity.⁶⁵ In the absence of experimental data specific to our system that could justify a more elaborate model, we adopt a simple and widely used pseudo-first-order kinetic scheme^65–67

where N_ads,eq(η) is the equilibrium number of adsorbed counterions on a single colloid and N_ads(t) is the number of adsorbed counterions at time t. For simplicity, we assume that the adsorption rate k is constant, independent of counterion concentration, and equal to the desorption rate. Assuming all counterions are monovalent, we can identify the number of adsorbed counterions with the adsorbed charge, N_ads,eq(η) = Z_ads,eq(η) and N_ads(t) = Z_ads(t). Following Diehl and Levin, we define the renormalized charge as Z* = Z − Z_ads,⁵⁵ where Z is the bare colloid charge arising from surface dissociation reactions, and −Z_ads represents the reduction due to counterion adsorption. In our system, the bare charge can be treated as constant, since the surface groups are strongly acidic and thus largely insensitive to local density variations, as shown in Fig. 1(a). Hence, we set dZ/dt = 0, or equivalently, Z = Z_eq. Substituting this relation into the kinetic equation yields a differential equation for the time evolution of the renormalized chargewhere τ = 1/k represents the characteristic timescale over which the colloidal surface charge relaxes towards equilibrium with its local chemical environment. Hence, the time evolution of Z* follows from first-order adsorption kinetics. To maintain consistency, we impose that the time evolution of the inverse screening length , as defined in eqn (7), mirrors that of the renormalized charge by making use of the one-to-one mapping between and , as is discussed in Section 3.1. Here, we follow a framework similar to that of Boon et al.,⁶⁸ who observed that the equilibration time of the electric double layer is typically much shorter than both the colloidal diffusion time and the timescales of most surface-charge reaction processes. This slow-reaction assumption allows us to describe the double layer in quasi-equilibrium. Consequently, we can employ the equilibrium DLVO potential at each instant, with the non-equilibrium effective charge Z*(t) and inverse screening length serving as time-dependent inputs. Furthermore, we note that the same differential equation eqn (12) has previously been used to model charging dynamics in aqueous nanochannels.⁶⁹

Fig. 1 Effect of charge regulation and charge renormalization on colloidal interactions. (a) Bare charge Z (dark blue) and renormalized charge Z* (light blue) as functions of packing fraction η. (b) Renormalized screening parameter as a function of η. (c) Total interaction potential βU(r) between two colloids as a function of center-of-mass distance r for two packing fractions, η = 0.0006 (blue) and 0.145 (red), representative of the initial gas and crystal phase, respectively. (d) Initial configuration used in the simulations consisting of a cubic crystallite surrounded by a dilute fluid phase. (e) Slice of thickness 5σ through the center of the simulation box, showing the renormalized charge Z* of individual particles (see colorbar for color coding). (f) Same slice as in (e), but with particle radii scaled by the local screening length according to 5/κσ + 0.7. Panels (e) and (f) correspond to time t = 2τ_d. The parameters used in the calculations are adopted from ref. 16.

Finally, we emphasize that our adsorption model assumes a uniform surface charge density on the colloids. This simplification contrasts with more detailed charge-regulating Poisson–Boltzmann calculations, which can yield heterogeneous and asymmetric charge distributions between colloids, and in some cases predict like-charge attractions.^70,71

2.3.3 Relaxation time τ for charge renormalization. Lastly, we provide a tentative estimate of the timescale τ for charge renormalization, acknowledging the inherent difficulties and therefore restricting ourselves to a lower bound. We focus on the relaxation of the innermost region of the electric double layer during a Brownian encounter—when two colloids approach closely enough for their double layers to overlap and become perturbed.⁷² Previous studies have shown that the diffuse part of the double layer relaxes rapidly and can be considered in near-instantaneous equilibrium under such conditions.^59,72,73 By contrast, the relaxation of condensed counterions involves a range of processes with different timescales and a strong dependence of interaction geometry.^58,72–75

In our simple adsorption model, the characteristic timescale τ in eqn (12) should be comparable to, or longer than, that of lateral counterion diffusion, in contrast to diffusion-limited models.^73,75 As a lower-bound estimate, we take τ to be the characteristic diffusion time along the colloid surface, obtained by dimensional analysis^72,73

where D_ion,c is the diffusion coefficient of ions in the condensed layer and a the colloid radius. It is worth emphasizing that the particular nature of the counterion association with the surface may result in τ ≫ τ_ion.

Using the diffusion coefficient of H⁺ in bulk water as an estimate for D_ion,c, we obtain a lower bound of τ_ion ≈ 0.01 ms. Experimental data suggests that the diffusion coefficient in the condensed layer may be comparable to its bulk value.^74,76,77 However, the structure of the electric double layer is complex, and the diffusion coefficients in the innermost layer remain debated.⁷⁸ In the presence of significant energy barriers attributed to lateral diffusion, the effective diffusion constant may be orders of magnitude smaller than in bulk.⁷³

Given these uncertainties, we also explore relaxation times τ that are significantly larger than τ_ion ≈ 0.01 ms. Additionally, systems with bulky counterions in nonpolar solvents, which are also actively studied,^24,54,79 may exhibit even slower relaxation dynamics, justifying our consideration of a broad range of τ values.

2.3.4 Generalization to charge regulation. The time-evolution equation for charge renormalization (eqn (12)) can also be interpreted as a simplified model for charge-regulating systems governed by ion adsorption,⁸⁰ where the rate-limiting timescale is determined by the kinetics of charge regulation rather than renormalization. In this case, the slow response is attributed to charge regulation dynamics, while the time dependence of charge renormalization itself is neglected. Charge regulation can be much slower than particle diffusion, sometimes requiring timescales on the order of minutes to hours to reach equilibrium.^{26–28,81,82} To reflect this broader applicability, we also consider relaxation times up to τ ≈ 3 s.

Here, we assume that the relaxation time τ = 1/k in the charging dynamics of eqn (12) is constant. In reality, the relaxation time τ is expected to increase with higher counterion concentrations [H⁺].⁸³ Nevertheless, the order of magnitude of τ is primarily determined by the specific system under consideration, as discussed in Section 2.3. For this study, we adopt a constant-τ assumption, which provides a reasonable approximation for the essential relaxation dynamics. We have added a discussion on the concentration dependence of τ and the rationale for employing the constant-τ model in Appendix C. Still, incorporating the concentration dependence of τ represents a natural extension for future work. Another natural extension would entail the decoupling of the time evolution of and Z*, as these can progress through different timescales.

3 Methods

3.1 Brownian dynamics with a dynamic colloid charge and screening length

In this work, we investigate the melting behavior of charged colloids using BD simulations. The time evolution of particle positions is evaluated using the Langevin integration scheme, where we follow the integration method proposed by Ermak and Yeh.⁸⁴ In this approach, the position r_i of particle i at time t + dt is given bywhere dt is the integration time step, and W(t) denotes a Wiener process that introduces stochastic noise corresponding to thermal fluctuations to the particles. The diffusion constant D determines the rate at which colloidal particles diffuse, setting the characteristic timescale τ_d = a²/D, known as the diffusion time τ_d. This timescale reflects how long it takes a particle to diffuse a distance comparable to its own size a. Finally, the interparticle forces are derived from the interaction potential U(r), which consists of a Yukawa potential accounting for the screened electrostatic repulsions and a hard-core WCA potential describing the excluded volume, as defined in eqn (1) and (2), respectively. By employing this integration scheme, we assume particles are in the overdamped limit and neglect hydrodynamic interactions. These approximations allow us to study macroscopic length and timescales. Nonetheless, hydrodynamic coupling under confinement¹⁸ and in the presence of charge regulation⁸³ is highly non-trivial and represents an interesting direction for extending our model.

In the simulations, we used a time step of dt = 8 × 10⁻⁶τ_d. The interaction potential was truncated and shifted to zero at a cutoff distance R_c = 8σ, following the truncated and shifted force method,⁸⁵ resulting in a negligible interaction strength of at most βU(R_c) ≈ 0.0001.

The interaction potential depends on the renormalized charge Z* and the inverse screening length ; we explained how to derive both of these parameters in Section 2. Here, the mean electrostatic potential ϕ_m for each packing fraction η was calculated using the Simpson integration method, while the PB cell model was solved numerically using the SCIPY SOLVE_BVP routine. Solutions for Z* and were obtained for the packing fraction range η ∈ [0, 0.735], and discretised into 4000 evenly spaced points. During the BD simulations, the values of Z* and were determined as continuous functions of η using linear interpolation.

To determine the instantaneous renormalized charge Z* and inverse screening length of each colloid during the simulation, we first compute its local packing fraction η using the Voronoi tesselation algorithm provided by voro++,⁴⁶ which is computationally demanding. Initially, the renormalized charge Z* and inverse screening length of each particle are set to their equilibrium values, and , for the initial configuration. To balance accuracy and efficiency, we update the equilibrium charge of each colloid via the Voronoi tesselation at a frequency of at least 1000 times per relaxation time τ. This ensures sufficient temporal resolution for systems with fast dynamics while reducing computational costs in more slowly relaxing systems. The charge Z*(t) and inverse screening length are updated at each simulation time step dt by integrating eqn (12). While a time-evolution equation for could in principle be derived from eqn (7), this would require knowledge of the time dependence of additional parameters in eqn (7) and would depend on the point of linearization. Instead, we evolve only Z*(t) and ensure consistency by using a one-to-one mapping between Z*(t) and at each time step. This mapping is valid within a specific range of packing fractions and allows us to accelerate the simulations. In the regime where a one-to-one mapping between Z* and exists, the value of consistent with the time-evolved Z* is determined via an effective packing fraction using linear interpolation. Simulations were carried out for a maximum duration of 160τ_d.

In eqn (1), a single renormalised inverse screening length is used, even though two interacting particles generally experience different local densities and thus different values. To account for this asymmetry, we approximate the interaction between particles i and j using the average inverse screening length, defined as . This approximation is expected to be reasonable in regions with small density gradients, such as those considered in Section 4.3. More detailed approaches could incorporate the full ion configurations around each colloid to define an effective screening length, but such treatments are beyond the scope of the present work.

3.2 Identification of crystalline particles

To determine when the colloidal crystals have fully melted, we use a machine-learning (ML) scheme to identify which particles are in a crystalline structure. To this end, we first classify the local environment of each particle by calculating the bond order parameters (BOPs).⁸⁶ These BOPs are based on spherical harmonics Y^m_l of order l, where m is an integer between m = −l and m = l. For each particle i, we calculatewhere r_ij = r_j − r_i is the distance vector between particle i and j, and the summation runs over all N_b(i) nearest neighbours of i. To distinguish more accurately between distinct thermodynamic phases, we calculate the averaged BOPS defined by Lechner and Dellago,⁸⁷ where we average q_lm(i) over the nearest neighbours and the particle itself⁸⁷

By averaging over −l ≤ m ≤ l, we define a rotationally invariant BOP

We calculated these BOPs (eqn (17)) for l ∈ [2, 11] using the PYTHON package PYSCAL3,⁸⁸ where the nearest neighbours were determined using a Voronoi tessellation.

The BOPs for all particles over the full simulation trajectory were projected onto two principal components using principal component analysis (PCA). A Gaussian mixture model (GMM) with two components was then fitted to this reduced dataset, under the assumption that one component represents the fluid phase and the other the FCC crystal phase. From the GMM, we obtained a crystallinity probability for each particle, which was summed across all particles and plotted as a function of time. To mitigate the influence of residual noise, we manually set the long-time crystallinity fraction to zero, effectively removing a small percentage of false positives.

In Appendix E, we provide a more detailed description of this method and demonstrate its robustness across various ML training conditions. The aim of this procedure is not to accurately identify the local structure of individual colloidal particles, but rather to capture the overall structural evolution of the system. Notably, these systems contain a substantial number of interface particles, which are notoriously difficult to classify.⁸⁹

4 Results and discussion

4.1 Effect of charge regulation and charge renormalization on colloidal interactions

Using the Poisson–Boltzmann cell theory as outlined in Section 2.2, we calculate both the bare colloidal charge Z, which arises from charge regulation, and the renormalized charge Z*, which accounts for charge renormalization, using eqn (6) and (7), respectively. Additionally, we calculate the corresponding renormalized inverse screening length using eqn (8). These calculations are based on a simple dissociation boundary condition with M = 10⁵ the number of surface groups, a = 465.5λ_B the colloid radius, κa = 1.185 the reduced screening parameter, and [H⁺]/K_a = 0.05 the ratio between the bulk concentration of hydrogen and the acid dissociation constant. All parameters except [H⁺]/K_a are taken from the 1997 paper by Larsen and Grier,¹⁶ where we use a Bjerrum length of λ_B = 0.7 nm. In Fig. 1(a), we plot the bare colloid charge Z (dark blue) and the renormalized charge Z* (light blue) as functions of the packing fraction η. We find that the bare charge Z remains approximately constant over a broad range of packing fractions but decreases sharply as η approaches the maximum packing fraction η_max. In contrast, the renormalized charge Z* exhibits a non-monotonic trend due to counterion condensation: it increases with η at low packing fractions and begins to decrease as Z drops near η_max. The increase in Z* with η has also been found in cell model calculations performed in the canonical ensemble.²⁴

Fig. 1(b) shows the corresponding renormalized inverse screening length . We clearly observe that increases rapidly with colloid packing fraction η in this low-salt system. This increase in results from the higher salt concentration at higher colloid density, as is described by the Gibbs-Donnan effect.⁹⁰ Our choice for the ratio between the bulk concentration of hydrogen and the acid dissociation constant [H⁺]/K_a was limited by the numerical stability of our calculations. However, in Appendix B we show that the solution to the PB equation in this parameter space regime is largely insensitive to the exact value of the bare surface charge. Therefore, the chosen value of [H⁺]/K_a is still expected to yield results that are representative of our model system.

The renormalized charge Z* and the inverse screening length can be employed in eqn (1) to define a dynamically evolving, density-dependent interaction potential. In Fig. 1(c), we illustrate this interaction at two representative packing fractions, η = 0.0006 and η = 0.145. This plot clearly shows that the Yukawa interaction is significantly more repulsive at low packing fraction. This reduction in repulsion at higher η is primarily due to enhanced electrostatic screening, reflected in the increase of . Although the renormalized charge Z* also rises with density, its effect is outweighed by the stronger screening.

4.2 Melting behavior of a cubic crystallite

Using the dynamically evolving, density-dependent Yukawa potential described above, we study the lifetime of a small cubic crystallite with an FCC structure, surrounded by a dilute fluid phase, for various values of the relaxation time τ. The initial configuration, shown in Fig. 1(d), consists of a cubic crystallite of 8788 particles arranged in an FCC structure with an initial packing fraction η_i, centered at the origin of the simulation box. This crystal is surrounded by 2000 randomly placed particles in a fluid phase. To achieve the desired global packing fraction η_g, we adjust the total volume of the simulation box while keeping the initial packing fraction η_i of the cubic crystal fixed. In all our simulations we use periodic boundary conditions in all three spatial directions. The parameters for all configurations used in this work, including η_g, η_i, and the relaxation time τ, are listed in Table 1, each identified by a code B_i with i ∈ [1, 16]. We first consider simulations B₉–B₁₂, in which the crystal has an initial packing fraction of η_i ≃ 0.231225, the global packing fraction is η_g ≃ 0.004189, and the fluid phase initially has a packing fraction of η ≃ 0.0015. Note that, at equilibrium, the colloids form a homogeneous fluid phase at a global packing fraction of η_g ≃ 0.004189.

Table 1 Parameters for the initial configurations of the BD simulations on a cubic crystallite with an FCC structure surrounded by a dilute fluid phase. The configurations are characterised by a code B_i, a global packing fraction η_g, the initial packing fraction of the crystal η_i, and the relaxation time τ used in the simulations

Code	η g	η i	τ/τd
B1	0.001047	0.057806	4
B2	0.001047	0.057806	0.4
B3	0.001047	0.057806	0.04
B4	0.001047	0.057806	0.004
B5	0.002094	0.115613	4
B6	0.002094	0.115613	0.4
B7	0.002094	0.115613	0.04
B8	0.002094	0.115613	0.004
B9	0.004189	0.231225	4
B10	0.004189	0.231225	0.4
B11	0.004189	0.231225	0.04
B12	0.004189	0.231225	0.004
B13	0.005236	0.404644	4
B14	0.005236	0.404644	0.4
B15	0.005236	0.404644	0.04
B16	0.005236	0.404644	0.004

---

In the dense crystalline region, the effective charge of the particles is higher than in the surrounding dilute fluid, as shown in Fig. 1(e). However, the particle interactions are less repulsive in the dense region due to the increased screening, reflected by a larger value of . This effect is illustrated in Fig. 1(f), where the colloid radius gives an impression of the interaction range. Interestingly, particles in the gas phase appear effectively larger to one another than those in the crystal, owing to the longer-range interactions at lower densities.

Typical configurations of the system at different stages of its melting process are shown in Fig. 2, illustrating how the expansion of the crystal depends on the relaxation time τ. At low values of τ, like τ = 0.04τ_d with τ_d the diffusion time of the colloids (Fig. 2(b)), the crystal rapidly expands into a dilute crystalline structure before eventually melting. In this regime, the charging mechanism responds almost instantaneously, so its time dependence is minimal. At longer relaxation times, τ = 4τ_d (Fig. 2(d)), the expansion proceeds much slower. After a time of t = 8τ_d, a dense FCC crystal still coexists with a surrounding dilute fluid phase. In contrast, for a short relaxation time of τ = 0.04τ_d, a large portion of the simulation box is already occupied by a dilute crystalline structure at the same time point.

Fig. 2 Melting behavior of a cubic crystal with an FCC structure. For clarity, only slices of thickness 6σ through the center of the simulation box are shown. The colors of the particles represent their renormalized charge Z*, as indicated by the color bar. (a) Slice of the initial configuration used in the melting simulations. (b) Typical configurations from the melting simulation with relaxation time τ = 0.04τ_d. (c) Time evolution of the fraction of particles identified as crystalline by our machine learning approach. (d) Typical configurations from the melting simulation with τ = 4τ_d.

Two main factors contribute to the stability of the crystal. The first arises from the relaxation dynamics within the dense region. Due to diffusion, particles tend to move on average from dense to dilute regions. Since the range of repulsion increases as the local density decreases, this can trigger a feedback loop that drives rapid crystal expansion. However, when the relaxation time τ is large, this feedback loop is suppressed because the particles have not yet adjusted to their new chemical environment. In other words, for small values of τ, the particles closely follow the light blue curves in Fig. 1(a) and (b). In contrast, for large τ, particles can effectively take a variety of paths through parameter space—spanning Z*λ_b/a, , and η—rather than remaining on an instantaneous response curve. This point is illustrated in Appendix D, where we study the time evolution of the charge and screening length of particles initially at the exterior and in the interior of the superheated crystal.

The second factor contributing to the stability of the crystal is the formation of an intermediate layer between the dense crystallite and the surrounding dilute fluid phase. When the charging dynamics is fast, i.e. small τ, particles at the surface of the crystallite quickly become more repulsive due to their lower local density. As a result, they experience strong outward repulsive forces that rapidly drive them away from the crystal center. In contrast, when the charging dynamics is slower, i.e. larger τ, surface particles do not immediately adjust to the reduced local density. Instead, they diffuse more gradually away from the crystallite, forming a dense, fluid-like layer. Although these particles have a lower local density than those inside the crystallite, their increased repulsion still plays a stabilizing role. This intermediate fluid layer acts as a barrier, exerting a counteracting force that resists further expansion of the crystal and thus enhances the stability and thereby the lifetime of the crystal. The strong repulsions in the intermediate fluid layer delay the breakup of the crystal, which ultimately melts due to weaker repulsions between colloids inside the crystal. Note that this effect strongly depends on the long screening length κ⁻¹ (resulting from the low ionic strength in the experiment), since the repulsions in the dilute regions need to be long-ranged in order to stabilise the crystal. The effects of increased stability of the crystal are shown in Fig. 2. At τ = 0.04τ_d, the crystal has completely dissolved by t = 40τ_d. In contrast, at τ = 4τ_d, crystalline regions remain visible at the same time point. The stabilizing fluid layer surrounding the crystal is also clearly visible in this figure.

To enable a more precise comparison between different values of τ, we quantify the extent of crystal melting by measuring the fraction of particles in an FCC-crystalline local environment. Fig. 2(c) shows the time evolution of the fraction of crystalline particles. The method used to identify the crystalline particles is described in Section 3.2 and detailed further in Appendix E. For each simulation, we calculated the BOPs for each particle. The datasets from simulations B₉–B₁₁ were projected onto the resulting principal components obtained by performing PCA on data from simulation B₁₂. Subsequently, a GMM was trained on the projected datasets. From Fig. 2(c), a clear distinction emerges between systems with a relaxation time τ > τ_d and those with τ < τ_d. For τ < τ_d, the crystalline fraction gradually decreases until there is no crystal left around t = 45τ_d. In contrast, for τ > τ_d, the fraction of crystalline particles initially plateaus, followed by a gradual decrease, with complete melting occuring at approximately t = 65τ_d.

In addition, we examine how the initial packing fraction η_i of the crystal influences its lifetime. We analyze three sets of simulations with varying initial packing fractions: B₁–B₄ with η_i = 0.057806, B₅–B₈ with η_i = 0.115613, and B₁₃–B₁₆ with η_i = 0.404644. For each set, PCA was performed using data from the simulation with τ = 0.004τ_d. The remaining three datasets in each set were projected onto the resulting principal components. A GMM was then trained and tested on each of the nine projected datasets. The resulting fraction of crystalline particles are presented in Fig. 3. Based on Fig. 3, we observe that increasing the initial packing fraction η_i enhances the lifetime of the crystal, at fixed relaxation time τ. More specifically, the lifetime of the crystal increases progressively from Fig. 3(a)–(c). This effect is particularly pronounced at larger relaxation times τ (see the red curves). As the initial packing fraction η_i increases, the difference between the two regimes τ < τ_d and τ > τ_d becomes increasingly evident. At low initial packing fraction η_i (Fig. 3(a)), the influence of τ on the melting rate is minimal. This is because density-dependent interaction effects are weak in dilute systems, and introducing a time lag in the response of the interactions has little impact. In contrast, at higher initial packing fractions η_i (Fig. 3(b) and (c)), the delayed response of the interactions plays a more important role. Here, larger values of τ noticeably extend the lifetime of the crystal. This suggests that as the system is driven further out of equilibrium, the crystal becomes more long-lived, and the delayed relaxation of the colloid surface chemistry plays an increasingly important role in the stabilization of the crystal.

Fig. 3 Melting behavior of cubic crystallites. Fraction of crystalline particles as a function of time t/τ_d for varying relaxation times τ and different initial packing fractions η_i of the crystal and global packing fraction η_g: (a) (η_i, η_g) = (0.057806, 0.001047) from simulations B₁, B₂, B₃, (b) (η_i, η_g) = (0.115613, 0.002094) from simulations B₅, B₆, B₇, and (c) (η_i, η_g) = (0.404644, 0.004189) from simulations B₁₃, B₁₄, B₁₅. In this dilute regime, the small variations in η_g between simulations have a negligible impact on the crystal lifetime.

To convert the dimensionless timescales discussed in this section into physical units, we refer to the particle diffusion coefficients reported by Larsen and Grier; either D = 0.12 μm² s⁻¹ (ref. 15) or D = 0.15 μm² s⁻¹.¹⁶ These values correspond to diffusion times of either τ_d = 0.89 s or τ_d = 0.71 s, respectively. Using the shorter diffusion time, τ_d = 0.71 s, the two relaxation times translate to τ = 0.04τ_d = 28 ms, and τ = 4τ_d = 2.8 s, and the lifetimes of the crystal range from 32 s–64 s.

In summary, for this specific set of initial conditions, the charge relaxation time—approximately ranging from milliseconds to seconds—strongly influences the lifetime of a superheated crystal, since a larger relaxation time can double the crystal lifetime. These results show that incorporating non-instantaneous charge regulation and charge renormalization effects into simulations of charged colloids enhances the stability and lifetime of crystalline structures.

4.3 Melting behavior of a crystal in a slab geometry

Next, we investigate a system with a geometry more closely resembling the experimental setup used by Larsen and Grier.^15,16 While the previous simulations of small compressed cubic crystals surrounded by a fluid phase allowed for a quantitative analysis of the effect of the relaxation time τ, the more complex geometry considered in the next simulations reveals that the specific details of the initial configuration also play a significant role in extending the lifetime of superheated crystals.

The initial configuration of the particles is shown in Fig. 4(a). An FCC crystal consisting of 22 000 particles with an initial packing fraction of η_i is positioned at the centre of an elongated simulation box. We apply periodic boundary conditions in all three spatial directions, resulting in two planar crystal-gas interfaces. The surrounding dilute gas phase contains 194 randomly placed particles. The total volume of the simulation box is adjusted to achieve a target global packing fraction η_g, while keeping η_i of the crystal fixed. The height of the box is maintained at approximately 31.5σ, corresponding to about 20.5 μm or 0.02 mm—matching the height of the experimental setup used by Larsen and Grier.¹⁵ All relevant parameters are provided in Table 2. A key difference between the simulations and experiments is, however, the use of periodic boundary conditions in the simulations. Furthermore, the simulations use a fully filled slab geometry, in contrast to the experiments, which consist of an unspecified number of layers. As a result, a direct quantitative comparison between the simulations and experiments is not feasible. Nevertheless, in the following, we highlight qualitative features that resemble those observed in the experimental system.

Fig. 4 Melting behavior of a crystal with an FCC structure in a slab geometry. (a) Initial configuration with particles colored according to their renormalized charge Z*, as shown by the color bar. The crystal has an initial packing fraction of η_i = 0.089215, and a global packing fraction of η_g = 0.009. The relaxation time is set to τ = 0.04τ_d. This corresponds to simulation P₃ of Table 2. (b) Time evolution of the system during expansion of the crystal. The particle colors represent the probability of belonging to the crystal phase. The images show slices of thickness 10σ, centered along the height of the simulation box with dimensions of approximately 72σ × 147σ. For visual clarity, all particle radii are set to 0.7σ.

Table 2 Parameters for the initial configurations of the BD simulations on a crystal in a slab geometry adjacent to a low-density gas phase, with the two phases separated by two planar interfaces. The configurations are characterised by a code P_i, a global packing fraction η_g, initial packing fraction of the crystal η_i, and the relaxation time τ used in the simulations

Code	η g	η i	τ/τd
P1	0.009	0.089215	4
P2	0.009	0.089215	0.4
P3	0.009	0.089215	0.04
P4	0.009	0.089215	0.004
P5	0.02	0.089213	4
P6	0.02	0.089213	0.4
P7	0.02	0.089213	0.04
P8	0.02	0.089213	0.004
P9	0.03	0.089213	4
P10	0.03	0.089213	0.4
P11	0.03	0.089213	0.04
P12	0.03	0.089213	0.004

---

We first examine a system with an initial packing fraction of the crystal of η_i = 0.089215, and a global packing fraction of η_g = 0.009. The relaxation time is set to τ = 0.04τ_d. However, given the low initial density of the superheated crystal, this value of τ is not expected to significantly influence the out-of-equilibrium dynamics (see Fig. 3). The simulation presented here corresponds to simulation P₃ in Table 2. To track the evolution of the crystal during expansion, we employ our machine learning method for identifying crystalline particles described in Section 3.2 and further detailed in Appendix E. To obtain the FCC class probabilities shown in Fig. 4, we projected the BOPs from P₃ onto the first two principal components obtained from simulation P₂. A GMM was then trained and tested on these reduced BOPs. Fig. 4(b) illustrates the expansion of the crystal, with particles color-coded according to their probability of belonging to the crystal phase. The color-coding was also validated against visual inspection of simulation trajectories in Supplementary Movies of the system.⁹¹ Due to the large system size, the simulation could not be run until it reaches full equilibrium.

We make the following observations when we run such a simulation. After an initial relaxation, the FCC crystal at the center of the box falls apart into irregularly shaped domains, separated by disordered particles, as shown in the left frame of Fig. 4(b). The subsequent frames illustrate the time evolution of these ordered regions. Notably, islands of FCC-structured particles persist and even grow in some cases in the central region of the box, where the local density remains relatively high. In these dense areas, ordered domains can span the full height of the simulation box. In contrast, the low-density regions near the left and right boundaries of the box, consist primarily of disordered particles. Between the dense crystal and the dilute phase, smaller crystalline domains are observed drifting through the disordered background. This drift is attributed to the slow, diffusion-driven expansion of the system, which gradually redistributes particles over time. Large ordered structures, similar to those shown in Fig. 4(b), persist throughout the entire duration of the simulations, which extended up to t = 135τ_d. These results indicate that the out-of-equilibrium coexistence between ordered and disordered phases described in this section has a macroscopic lifetime of at least approximately 1.5 min, and potentially significantly longer.

In addition, we also consider conditions analogous to those in the experimental study, where the system equilibrated into a disordered fluid with a global packing fraction of η_g = 0.02. However, in our simulations at this global packing fraction, the crystal did not melt. This behavior is shown in the Supplementary Video⁹¹ and in Fig. 5(a), which shows the expansion of the crystal from simulation P₇. The color-coding of the particles denote the FCC class probability, determined using the method described in Section 3.2, with principal components derived from simulation P₆.

Fig. 5 Expansion of a crystal in a slab geometry at a global packing fraction η_g = 0.02, corresponding to simulation P₇ of Table 2. (a) The first four frames show the evolution of the crystalline regions during expansion of the crystal. (b) Packing fraction gradient across the simulation box at time t = 64τ_d, with the color bar indicating the local packing fraction η.

From the final two frames of Fig. 5(a), it can be inferred that the particles have reached the boundaries of the simulation box. However, the system did not show any signs of melting throughout the entire simulation. This suggests that at a global packing fraction of η_g = 0.02, our simulations do not reproduce the disordered equilibrium fluid observed experimentally by Larsen and Grier.^15,16 This discrepancy suggests that the effective strength or range reported for the Yukawa potential in eqn (1) in the experimental system may be overestimated.

The first two frames of Fig. 5(a) show the system before the colloids reach the boundaries of the simulation box. During this early stage, particles located near the dilute region are more likely to be assigned to the FCC class probability than particles farther inward. This is an artifact of the classification method. Once the system has fully expanded and the particles form a continuous structure – due to periodic boundary conditions – this artifact vanishes. Fig. 5(b) displays the packing fraction gradient of simulation P₇ at t = 64τ_d. We observe that the only significant density variation occurs along the direction of expansion.

While the shapes of the ordered regions in Fig. 4(b) and 5(b) resemble those observed in the experimental system, we did not observe any significant density differences between the ordered and disordered regions. This suggest that additional factors are likely necessary to explain the coexistence of dense crystalline and dilute fluid regions reported by Larsen and Grier.^15,16 Nonetheless, our simulations display long-lived order–disorder coexistence, characterized by islands of FCC-ordered particles that evolve in shape and drift as a result of the expansion of the system. These results indicate that not all features observed in experiments necessarily stem from like-charge attraction. Moreover, our findings highlight the strong effect of system size and geometry on the stability and lifetime of superheated crystals. Although the compressed cubic crystallites discussed in Section 4.2 were initialized at nearly ten times the density of the crystal in the slab geometry, the latter maintained ordered regions for substantially longer durations.

4.4 Walls and like-charge attractions

One aspect not addressed in our analysis is the complex interplay between walls and colloids. To place our long-box simulations in the broader context of wall-induced like-charge attraction, we briefly discuss previous studies examining electrostatic and hydrodynamic effects on the effective interactions of charged colloids near walls.

Many experimental studies have investigated the effective interaction potential between pairs of colloids, positioned either near a single plate or confined between two glass plates.^16,92–103 Several of these studies report indications of like-charge attraction,^{16,92,94,95,97,98,100–103} although possible experimental artifacts have also been raised.^{92,93,100,104} The most compelling evidence for an attractive potential between polystyrene particles arises in samples confined between two glass walls,^{94,95,98,100–103} where strong confinement was found to be essential.^94,95,101 Studies based on Poisson–Boltzmann theory showed, however, that the interaction potential between two like-charged colloids confined within a cylinder remains strictly repulsive.^105–108 Importantly, both the strongly confined geometry and the cylindrical confinement differ significantly from that considered in the present work.^15,16 The effective interaction between colloids near a single wall has also been studied, where it was suggested that planes of colloids may generate partial confinement, allowing particles to act as walls for one another and transmit the influence of the container into the bulk crystal.^16,101 While Larsen and Grier¹⁶ reported an attractive interaction potential for polystyrene colloids near a single wall, these findings have been challenged by subsequent studies,^94,96,97 leaving the existence of such attractions for polystyrene colloids unclear. In contrast, attractive interactions have been more consistently observed for weakly acidic silica colloids.⁹² Interestingly, for silica systems, the influence of a single wall decreases with decreasing ionic strength,^92,95,99 whereas low-salt conditions are required for stabilizing metastable crystals.^15,16 Thus, for polystyrene colloids, confinement-induced like-charge attractions have been observed only under strong confinement, whereas for silica colloids, the measured wall-induced attractions decrease with decreasing salt concentration—contrasting with the low-salt conditions required to maintain the metastability of superheated crystals. These observations suggest that confinement-induced attractions cannot explain the metastability of superheated crystals formed by polystyrene latex colloids.

Finally, Squires and Brenner¹⁸ showed that non-equilibrium hydrodynamic coupling between a pair of colloids located at a distance h from a container wall can give rise to an apparent like-charge attraction. Their analysis assumed constant, uniform negative charges on all three objects, with the plate's surface charge density treated as a fitting parameter.

Notably, recent observations of like-charge attraction in samples of large membrane-coated colloids revealed an asymmetry with respect to the sign of the colloid charge: only negatively charged colloids exhibited particle clustering.^109,110 Further experimental work showed that negatively charged particles can display like-charge attraction in polar solvents, whereas positively charged colloids do so in nonpolar solvents.¹⁴ Kubincová et al. attributed this asymmetry to an electrosolvation effect,¹³ a conjecture that has since been supported by further studies reporting promising evidence.^111–113 Nonetheless, it remains an open question whether the cases of like-charge attraction discussed here arise from the same underlying mechanisms.

5 Summary and conclusions

In this study, we have performed simulations of charged colloids, in which the colloid charge was considered as a dynamically evolving, density-dependent parameter. Using Poisson–Boltzmann cell theory calculations, we determined the renormalized and regulated charge as functions of packing fraction. To capture the non-instantaneous charging behavior of colloidal surfaces, we introduced a time-dependent differential equation governing the charging dynamics during the simulations. This approach was then applied to low-salt suspensions of superheated colloidal crystals. To the best of our knowledge, this represents the first simulation implementation of charge-regulating colloids with explicitly density-dependent charges.

The first series of simulations was initialised by placing a cubic colloidal crystal next to a colloidal fluid. Although all simulated crystals eventually melted, their lifetimes were significantly extended when both density- and time-dependence were incorporated into the determination of colloid charge and screening length. More specifically, the simulations revealed two distinct regimes based on the relaxation time τ, relative to the diffusion timescale τ_d: τ > τ_d and τ < τ_d. Systems with τ > τ_d exhibited notably longer crystal lifetimes and are most relevant for experimental settings where charge regulation is slow. In contrast, systems with τ < τ_d, where the relaxation dynamics are effectively instantaneous, displayed faster melting dynamics. Within the context of our experimental model system, we observed an out-of-equilibrium crystal-fluid coexistence lasting at least 46 s. Additionally, increasing the initial packing fraction of the crystal further enhanced its stability and lifetime.

In the second set of simulations, a colloidal crystal slab was positioned adjacent to a low-density colloidal gas phase, with the two phases separated by two planar interfaces as opposed to six interfaces in the case of a cubic crystal. As the crystal expanded, a dynamic coexistence emerged between crystalline and fluid-like regions, with ordered crystalline domains drifting through surrounding disordered regions of colloids. These crystalline regions grew in some areas and dissolved in others, leading to a continuously evolving structure. This non-equilibrium phase coexistence persisted throughout the entire simulation duration of t = 135τ_d, closely resembling the heterogeneous structures observed by Larsen and Grier.¹⁶ However, in contrast to the experimental observations, our simulations did not exhibit significant density differences between the crystalline and fluid regions, as both phases maintained comparable local densities. This discrepancy suggests that additional physical mechanisms are necessary to fully account for the observed metastability and phase coexistence in charged colloidal crystals. Furthermore, our results highlight the significant influence of system size and geometry on the stability and lifetime of superheated crystals. Although the compressed cubic crystallites discussed in Section 4.2 were initialized at nearly ten times the density of the crystal in the slab geometry, the latter sustained ordered regions for substantially longer timescales.

Altogether, incorporating dynamically evolving, density-dependent colloidal interactions leads to unexpected dynamic melting behavior, characterized by long-lived coexistence between crystalline domains and disordered regions. These results underscore the importance of conducting more detailed simulations of systems with density-dependent charging, as such interactions can give rise to non-intuitive and emergent out-of-equilibrium phenomena. Notably, our findings depend on the initial colloid configuration and the charge-regulation curve used. Changing the surface-charging mechanism under scrutiny is therefore expected to produce markedly different out-of-equilibrium responses, offering intriguing possibilities for tuning colloidal phase behaviour. Finally, we note that metastable crystalline clusters are often explained in terms of effective like-charge attractions. In contrast, our results show that the lifetime of metastable crystals can be significantly extended by the slow relaxation dynamics of surface charges—arising from charge regulation and charge renormalization—while the colloidal interactions remain purely repulsive. It would be interesting to further investigate the kinetic relaxation mechanism presented here in combination with more detailed charge-regulating Poisson–Boltzmann calculations, in order to capture the influence of solvent effects and heterogeneous surface charge distributions on like-charged colloidal interactions.

Author contributions

Laura Jansen: conceptualization (equal); data curation (equal); formal analysis (lead); investigation (lead); methodology (equal); software (lead); visualization (lead); writing – original draft (equal); writing – review & editing (supporting). Thijs ter Rele: conceptualization (equal); data curation (equal); investigation (supporting); methodology (equal); software (supporting); visualization (supporting); writing – original draft (equal); writing – review & editing (supporting). Marjolein Dijkstra: conceptualization (equal); funding acquisition (lead); methodology (equal); project administration (lead); supervision (lead); writing – review & editing (lead).

Conflicts of interest

There are no conflicts to declare.

Data availability

The code used in generating and analyzing the data presented in work are included in the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5sm00713e.

The data itself is not included, following from the large size of the generated files. However, it is readily available on request.

Acknowledgements

M. D. and T. t. R acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. ERC-2019-ADG 884902 SoftML). We thank René van Roij for useful discussions.

Notes and references

1. E. Matijevic, Medical applications of colloids, Springer Science & Business Media, 2008.
2. E. J. W. Verwey and J. T. G. Overbeek, Theory of the stability of lyophobic colloids: the interaction of sol particles having and electric double layer, Elsevier Publishing Company, 1948.
3. B. Derjaguin and L. Landau, Prog. Surf. Sci., 1993, 43, 30–59.
4. J. N. Israelachvili, Intermolecular and Surface Forces, Elsevier, 3rd edn, 2011.
5. G. Trefalt, T. Palberg and M. Borkovec, Curr. Opin. Colloid Interface Sci., 2017, 27, 9–17.
6. T. Yoshiyama, I. Sogami and N. Ise, Phys. Rev. Lett., 1984, 53, 2153–2156.
7. H. Yoshida, K. Ito and N. Ise, J. Am. Chem. Soc., 1990, 112, 592–596.
8. K. Ito, H. Nakamura, H. Yoshida and N. Ise, J. Am. Chem. Soc., 1988, 110, 6955–6963.
9. B. V. R. Tata, M. Rajalakshmi and A. K. Arora, Phys. Rev. Lett., 1992, 69, 3778–3781.
10. K. Ito, H. Yoshida and N. Ise, Science, 1994, 263, 66–68.
11. H. Yoshida, N. Ise and T. Hashimoto, J. Chem. Phys., 1995, 103, 10146–10151.
12. B. V. R. Tata, E. Yamahara, P. V. Rajamani and N. Ise, Phys. Rev. Lett., 1997, 78, 2660–2663.
13. A. Kubincová, P. H. Hünenberger and M. Krishnan, J. Chem. Phys., 2020, 152, 104713.
14. S. Wang, R. Walker-Gibbons, B. Watkins, M. Flynn and M. Krishnan, Nat. Nanotechnol., 2024, 19, 485–493.
15. A. E. Larsen and D. G. Grier, Phys. Rev. Lett., 1996, 76, 3862–3865.
16. A. E. Larsen and D. G. Grier, Nature, 1997, 385, 230–233.
17. T. Palberg and M. Würth, Phys. Rev. Lett., 1994, 72, 786.
18. T. M. Squires and M. P. Brenner, Phys. Rev. Lett., 2000, 85, 4976–4979.
19. C. P. Royall, M. E. Leunissen, A.-P. Hynninen, M. Dijkstra and A. Van Blaaderen, J. Chem. Phys., 2006, 124, 244706.
20. T. Kanai, N. Boon, P. J. Lu, E. Sloutskin, A. B. Schofield, F. Smallenburg, R. Van Roij, M. Dijkstra and D. A. Weitz, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2015, 91, 030301(R).
21. T. Markovich, D. Andelman and R. Podgornik, Europhys. Lett., 2016, 113, 26004.
22. B. W. Ninham and V. A. Parsegian, J. Theor. Biol., 1971, 31, 405–428.
23. J. C. Everts, N. Boon and R. Van Roij, Phys. Chem. Chem. Phys., 2016, 18, 5211–5218.
24. J. E. Hallett, D. A. J. Gillespie, R. M. Richardson and P. Bartlett, Soft Matter, 2018, 14, 331–343.
25. B. Zoetekouw, Phase behavior of charged colloids: many-body effects, charge renormalization and charge regulation, 2006, https://dspace.library.uu.nl/handle/1874/12536.
26. S. Biswas and D. Chattoraj, J. Colloid Interface Sci., 1998, 205, 12–20.
27. M. J. Avena and L. K. Koopal, Environ. Sci. Technol., 1999, 33, 2739–2744.
28. L. Krumina, J. P. Kenney, J. S. Loring and P. Persson, Chem. Geol., 2016, 427, 54–64.
29. G. Frens and J. Overbeek, J. Colloid Interface Sci., 1972, 38, 376–387.
30. S. Alexander, P. M. Chaikin, P. M. Grant, G. J. Morales, P. Pincus and D. Hone, J. Chem. Phys., 1984, 80, 5776–5781.
31. E. Trizac, L. Bocquet, M. Aubouy and H. H. Von Grünberg, Langmuir, 2003, 19, 4027–4033.
32. L. Belloni, Colloids Surf., A, 1998, 140, 227–243.
33. T. Gisler, S. Schulz, M. Borkovec, H. Sticher, P. Schurtenberger, B. D’Aguanno and R. Klein, J. Chem. Phys., 1994, 101, 9924–9936.
34. M. Quesada-Pérez, J. Callejas-Fernández and R. Hidalgo-Álvarez, Colloids Surf., A, 1999, 159, 239–252.
35. T. Palberg, W. Mönch, F. Bitzer, R. Piazza and T. Bellini, Phys. Rev. Lett., 1995, 74, 4555–4558.
36. H. Von Grünberg, J. Colloid Interface Sci., 1999, 219, 339–344.
37. L. F. Rojas-Ochoa, R. Castañeda-Priego, V. Lobaskin, A. Stradner, F. Scheffold and P. Schurtenberger, Phys. Rev. Lett., 2008, 100, 178304–178307.
38. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
39. L. Filion, R. Ni, D. Frenkel and M. Dijkstra, J. Chem. Phys., 2011, 134, 134901.
40. M. Deserno and C. Holm, Cell Model and Poisson-Boltzmann Theory: A Brief Introduction, Springer, Dordrecht, 2001, pp. 27–52.
41. F. Smallenburg, N. Boon, M. Kater, M. Dijkstra and R. Van Roij, J. Chem. Phys., 2011, 134, 074505.
42. H. Ohshima, Electrical phenomena at interfaces and biointerfaces, John Wiley & Sons, 2012.
43. D. J. Griffiths, Introduction to Electrodynamics, Cambridge University Press, 2017.
44. P. Biesheuvel and M. A. C. Stuart, Langmuir, 2004, 20, 2785–2791.
45. M. E. Brito, G. Nägele and A. R. Denton, J. Chem. Phys., 2023, 159, 204904.
46. C. H. Rycroft, Chaos, 2009, 19, 041111.
47. G. S. Manning, J. Chem. Phys., 1969, 51, 924–933.
48. G. S. Manning, Acc. Chem. Res., 1979, 12, 443–449.
49. D. Hinderberger, H. W. Spiess and G. Jeschke, Appl. Magn. Reson., 2009, 37, 657–683.
50. T. S. Lo, B. Khusid and J. Koplik, Phys. Rev. Lett., 2008, 100, 128301.
51. J. Granot, J. Feigon and D. R. Kearns, Biopolymers, 1982, 21, 181–201.
52. G. S. Manning, J. Phys. Chem. B, 2007, 111, 8554–8559.
53. Q. Tang and M. Rubinstein, Soft Matter, 2022, 18, 1154–1173.
54. D. A. J. Gillespie, J. E. Hallett, O. Elujoba, A. F. C. Hamzah, R. M. Richardson and P. Bartlett, Soft Matter, 2014, 10, 566–577.
55. A. Diehl and Y. Levin, J. Chem. Phys., 2004, 121, 12100–12103.
56. V. Sanghiran and K. S. Schmitz, Langmuir, 2000, 16, 7566–7574.
57. S. Shulepov, S. Dukhin and J. Lyklema, J. Colloid Interface Sci., 1995, 171, 340–350.
58. P. M. Biesheuvel, Langmuir, 2002, 18, 5566–5571.
59. S. S. Dukhin and J. Lyklema, Faraday Discuss. Chem. Soc., 1990, 90, 261.
60. M. Muthukumar, J. Chem. Phys., 2004, 120, 9343–9350.
61. S. Ghosh and A. Kundagrami, J. Chem. Phys., 2024, 160, 084909.
62. Y. Shi, H. Peng, J. Yang and J. Zhao, Macromolecules, 2021, 54, 4926–4933.
63. K. Radhakrishnan and S. P. Singh, Phys. Rev. E, 2022, 105, 044501.
64. R. Nikam, X. Xu, M. Kanduč and J. Dzubiella, J. Chem. Phys., 2020, 153, 044904.
65. Y. S. Ho, J. C. Ng and G. McKay, Sep. Purif. Methods, 2000, 29, 189–232.
66. J.-P. Simonin, Chem. Eng. J., 2016, 300, 254–263.
67. W. Rudzinski and W. Plazinski, J. Phys. Chem. B, 2006, 110, 16514–16525.
68. W. Boon, M. Dijkstra and R. van Roij, Phys. Rev. Lett., 2023, 130, 058001.
69. T. M. Kamsma, W. Q. Boon, T. ter Rele, C. Spitoni and R. van Roij, Phys. Rev. Lett., 2023, 130, 268401.
70. A. Majee, M. Bier, R. Blossey and R. Podgornik, Phys. Rev. Res., 2020, 2, 043417.
71. H. Ruixuan, A. Majee, J. Dobnikar and R. Podgornik, Eur. Phys. J. E:Soft Matter Biol. Phys., 2023, 46, 115.
72. J. Lyklema, Pure Appl. Chem., 1980, 52, 1221–1227.
73. D. Weaver and D. Feke, J. Colloid Interface Sci., 1985, 103, 267–283.
74. J. Lyklema, H. Van Leeuwen and M. Minor, Adv. Colloid Interface Sci., 1999, 83, 33–69.
75. J. W. Krozel, J. Colloid Interface Sci., 1994, 163, 437–453.
76. J. Lyklema, Mol. Phys., 2002, 100, 3177–3185.
77. J. Lyklema, S. Rovillard and J. De Coninck, Langmuir, 1998, 14, 5659–5663.
78. R. Hartkamp, A. Biance, L. Fu, J. Dufrêche, O. Bonhomme and L. Joly, Curr. Opin. Colloid Interface Sci., 2018, 37, 101–114.
79. G. Hussain, A. Robinson and P. Bartlett, Langmuir, 2013, 29, 4204–4213.
80. G. Trefalt, S. H. Behrens and M. Borkovec, Langmuir, 2015, 32, 380–400.
81. D. Lis, E. H. G. Backus, J. Hunger, S. H. Parekh and M. Bonn, Science, 2014, 344, 1138–1142.
82. B. L. Werkhoven, J. C. Everts, S. Samin and R. Van Roij, Phys. Rev. Lett., 2018, 120, 264502.
83. K. Takae and H. Tanaka, Soft Matter, 2018, 14, 4711–4720.
84. D. Ermak and Y. Yeh, Chem. Phys. Lett., 1974, 24, 243–248.
85. M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford University Press, 2nd edn, 2017.
86. P. J. Steinhardt, D. R. Nelson and M. Ronchetti, Phys. Rev. B:Condens. Matter Mater. Phys., 1983, 28, 784–805.
87. W. Lechner and C. Dellago, J. Chem. Phys., 2008, 129, 114707.
88. S. Menon, G. Leines and J. Rogal, J. Open Source Software, 2019, 4, 1824.
89. W. Gispen, A. P. De Alba Ortíz and M. Dijkstra, J. Chem. Phys., 2025, 162, 024502.
90. F. G. Donnan, Z. Elektrochem. Angew. Phys. Chem., 1911, 17, 572–581.
91. SI for the Paper, https://pubs.rsc.org/en/Content/ArticleLanding/2025/SM/D5SM00713E.
92. M. Polin, D. G. Grier and Y. Han, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 041406.
93. J. Baumgartl, J. L. Arauz-Lara and C. Bechinger, Soft Matter, 2006, 2, 631.
94. D. G. Grier and Y. Han, J. Phys.: Condens. Matter, 2004, 16, S4145–S4157.
95. Y. Han and D. G. Grier, Phys. Rev. Lett., 2003, 91, 038302.
96. C. Franck, M. Covelli and R. V. Durand, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2003, 67, 041402.
97. M. Brunner, C. Bechinger, W. Strepp, V. Lobaskin and H. H. Von Grünberg, Europhys. Lett., 2002, 58, 926–965.
98. A. Ramírez-Saito, M. Chávez-Páez, J. Santana-Solano and J. L. Arauz-Lara, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2003, 67, 050403.
99. S. H. Behrens and D. G. Grier, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2001, 64, 050401.
100. K. S. Rao and R. Rajagopalan, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 3227–3233.
101. J. C. Crocker and D. G. Grier, Phys. Rev. Lett., 1996, 77, 1897–1900.
102. M. D. Carbajal-Tinoco, F. Castro-Román and J. L. Arauz-Lara, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 53, 3745–3749.
103. G. M. Kepler and S. Fraden, Phys. Rev. Lett., 1994, 73, 356–359.
104. M. Gyger, F. Rückerl, J. Käs and J. Ruiz-García, J. Colloid Interface Sci., 2008, 326, 382–386.
105. J. E. Sader and D. Y. C. Chan, Langmuir, 1999, 16, 324–331.
106. J. E. Sader and D. Y. Chan, J. Colloid Interface Sci., 1999, 213, 268–269.
107. J. C. Neu, Phys. Rev. Lett., 1999, 82, 1072–1074.
108. E. Trizac and J.-L. Raimbault, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 6530–6533.
109. M. M. Baksh, M. Jaros and J. T. Groves, Nature, 2004, 427, 139–141.
110. E. W. Gomez, N. G. Clack, H.-J. Wu and J. T. Groves, Soft Matter, 2009, 5, 1931.
111. S. Wang, R. Walker-Gibbons, B. Watkins, B. Lin and M. Krishnan, Nat. Commun., 2025, 16, 2872.
112. R. Walker-Gibbons, A. Kubincová, P. H. Hünenberger and M. Krishnan, J. Phys. Chem. B, 2022, 126, 4697–4710.
113. A. Behjatian, R. Walker-Gibbons, A. A. Schekochihin and M. Krishnan, Langmuir, 2022, 38, 786–800.

---