Quasi-two-dimensional dispersions of Brownian particles with competitive interactions: phase behavior and structural properties

Abstract

Competing short-range attractive (SA) and long range repulsive (LR) particle interactions can be used to describe three-dimensional charge-stabilized colloid or protein dispersions at low added salt concentrations, as well as membrane proteins with interaction contributions mediated by lipid molecules. Using Langevin dynamics (LD) simulations, we determine the generalized phase diagram, cluster shapes and size distributions of a generic quasi-two-dimensional (Q2D) dispersion of spherical SALR particles confined to in-plane motion inside a bulk fluid. The SA and LR interaction parts are modelled by a generalized Lennard-Jones potential and a screened Coulomb potential, respectively. The microstructures of the detected equilibrium and non-equilibrium Q2D phases are distinctly different from those observed in three-dimensional (3D) SALR systems, by exhibiting different levels of hexagonal ordering. We discuss a thermodynamic perturbation theory prediction for the metastable binodal line of a reference system of particles with SA interactions only, which in the explored Q2D-SALR phase diagram region separates cluster from non-clustered phases. The transition from the high-temperature (small SA) dispersed fluid (DF) phase to the lower-temperature equilibrium cluster (EC) fluid phase is characterised by a low-wavenumber peak height of the static structure factor (corresponding to a thermal correlation length of about twice the particle diameter) featuring a distinctly smaller value (≈1.4) than in 3D SALR systems. With decreasing temperature (increasing SA), the cluster morphology changes from disk-like shapes in the equilibrium cluster phase, to double-stranded anisotropic hexagonal cluster segments formed in a cluster-percolated (CP) gel-like phase. This transition can be quantified by a hexagonal order parameter distribution function. The mean cluster size and coordination number of particles in the CP phase are insensitive to changes in the attraction strength.

---

1 Introduction

During the past two decades, strong attention has been paid to Brownian particle dispersions with competitive short-range attractive (SA) and long-range repulsive (LR) inter-particle forces.^1–3 The competition between attraction and repulsion makes these excellent model systems for exploring the physics of copious soft matter and biological systems, including charged-stabilized colloids with added depletants,^1,4–8 lysozyme protein solutions under low-salinity conditions,^1,9–14 Y-shaped monoclonal antibodies,^15–17 and the reentrant liquid condensation of ribonucleoprotein-RNA complexes,¹⁸ to name just a few. Adjusting the respective influence of SA and LR forces allows to control the microstructure of the dispersions, which is relevant to applications such as in the formulation of pharmaceutical drugs and protein crystallization.² Owing to their importance both for fundamental theory and in industrial applications, the structure and phase behavior of three-dimensional (bulk) dispersions of SALR particles have been extensively investigated theoretically, numerically and experimentally. For example, integral equation theory,^19,20 discrete perturbation theory (DPT),^21–23 density functional theory,^24–27 and computer simulations^23,28 have been used to study the phase behavior of 3D-SALR systems. In ref. 29 and the mini-review,³⁰ Monte Carlo simulations in three dimensions (3D) have been reported for a SALR potential consisting of hard-core and short-range attractive square-well parts, augmented by a long-range linear repulsive ramp potential. In ref. 31 a phase diagram in two dimensions (2D) is obtained for a triangular lattice model where neighbouring particles are attractive and next-neighbour particles are repulsive. Here, a cluster phase is found, as well as stripe-like phases. Recent 3D-MD simulations of SALR particles in between two parallel walls with one wall being attractive³² have shown flattened clusters or stripes, respectively, formed by the adsorbed particles. Generic SALR pair interaction potentials have been also employed to mimick the behavior of aqueous lysozyme solutions for low added salt concentrations.^13,14

In contrast to pre-dominantly attractive particles fluids where macroscopic gas–liquid type phase separation is observed, the repulsive interaction part in SALR suspensions primarily stints the size of growing clusters, thus frustrating macroscopic phase separation.^33–35 Due to the strongly differing length scales of attraction and repulsion, there is a rich phase behavior of 3D-SALR systems including equilibrium and gel-like cluster phases.^23,36–38 Using Monte-Carlo simulations for two generic SALR potentials, Godfrin et al.²³ have mapped out a generalized 3D phase diagram including four different phases distinguishable by their respective cluster size distribution (CFD) functions, namely equilibrium dispersed fluid (DF) and equilibrium cluster (EC) states, and random percolated (RP) and cluster percolated (CP) states. In the temperature versus particle concentration state diagram, the latter two states are separated from the former ones by a percolation line, while the EC and CP cluster states (non-clustered DF and RP states) are located below (above) the binodal line of a short-range attractive (SA) reference potential system. As argued by Ruiz-Franco and Zaccarrelli,³ the phase states observed by Godfrin et al. are located well above a so-called λ-line in the temperature-concentration state diagram representing systems for which the static structure factor diverges at a small but non-zero wavenumber q_λ > 0. The λ-line marks the onset of a low-temperature region where the system becomes unstable to periodic density fluctuations of long wavelength 2π/q_λ, giving rise to ordered modulated phases such as stripe or lamellar phases, self-organized from rather monodisperse clusters. We note that microphase separation into ordered modulated phases is not limited to colloidal systems, but also occurs in micellar systems where the micelles exhibit SALR-type interactions, as shown in ref. 39.

In contrast, the higher-temperature disordered EC and CP cluster phases identified by Godfrin et al. consist of finite-sized clusters with a broad size distribution. A rough estimate of the characteristic wavenumbers q_c and q_λ (λ-line) can be obtained using a simple random phase approximation calculation of the SALR static structure factor, with a hard-sphere or hard-disk reference potential in the 3D and 2D cases, respectively.²⁷

Regarding the (metastable) binodal of the SA reference potential system separating the EC and CP phases from the higher-temperature DF and RP phases, the Noro–Frenkel extended law of corresponding states (ELCS) asserts that the thermodynamics (and equilibrium structure) of all three-dimensional SA systems is practically the same, independent of the details of the underlying SA potential, provided its attraction range is less than 25% of the particle diameter, for systems of equal temperature (attraction strength), second virial coefficient and volume fraction, expressed in units of the according critical point values at the liquid–liquid phase separation.⁴⁰ For protein solutions where the electrostatic repulsion is strongly screened, it has been shown experimentally that near to the liquid–liquid phase separation, the ELCS behavior further extends to collective diffusion.⁴¹ While the ELCS does not hold as such for 3D-SALR systems,³ an universal state diagram is obtained, provided temperature and particle concentration are scaled by the critical point temperature and concentration of the SA reference system, respectively.²³

As discussed in ref. 2 and 10, the (onset of) finite-sized particle clusters formation is globally indicated by a prepeak of the static structure factor, S(q), at a wavenumber, q = q_c, distinctly smaller than the wavenumber, q_m, of the principal peak associated with the nearest neighbor shell of particles. In contrast to q_m, the prepeak position wavenumber q_c is in general rather insensitive to changes in the particles concentration. The prepeak of S(q) is also commonly referred to as intermediate range order peak, to emphasize that its existence does not necessarily indicate a genuine clustered fluid phase.^2,10 While the present work is not concerned with the dynamics of SALR systems, we note that a prepeak is theoretically predicted also for the hydrodynamic function, H(q), of 3D-SALR systems characterizing collective diffusion,^28,42 and this feature has been experimentally detected for lyzosyme solutions under low-salt conditions in the DF phase region.¹⁴

For 3D-SALR systems, a notable effort was made to identify simple criteria characterizing the transition from the DF to the EC equilibrium phase. From their MC simulations, Godfrin et al.²³ conclude that this transition is hallmarked by a prepeak height value S(q_c) ≈ 2.7. This value is accidentally close to the Hansen–Verlet freezing transition value S(q_m) ≈ 2.85 for the height of the static structure factor principal peak of a three-dimensional hard-sphere system, located at the wavenumber q_m^43,44 associated with the nearest neighbor shell of particles. The Hansen–Verlet freezing transition criterion can be mapped, incidentally, to a dynamic criterion in terms of the long-time self-diffusion coefficient, both for 3D and 2D systems.⁴⁵ The freezing transition value is enlarged slightly to S(q_m) = 3.1 for low-salt dispersions of highly charged colloids, in which case it signals a fluid to bcc crystal phase transition.⁴⁶

Observation of the prepeak height value S(q_c) = 2.7 is not sufficient to unambiguosly indicate a DF to EC phase transition in three dimensions. According to Bollinger and Truskett³³ (see also ref. 3 and 47), this condition should be amended by the condition 2.0 ≲ ξ_T/σ ≲ 3.0 for the thermal correlation length ξ_T in units of the particle diameter σ, where ξ_T is identified as the inverse half-width of a quadratic fit of the inverse structure factor prepeak around q_c. As we will discuss also in the context of Q2D dispersions, fluid-like systems with a distinct cluster peak S(q_c) exhibit a damped, long-distance oscillatory part of the radial distribution function g(r) of characteristic wavelength 2π/q_c. Additionally to the prepeak-associated indicators S(q_c) and ξ_T, an empirical relation between particle coordination number and mean cluster size has been identified as being likewise indicative of the DF to EC transition.²³ Since in the present study only the depth of the short-ranged attractive part of the potential is varied, phases additional to the ones reported here may appear, e.g., for stronger attraction strengths and a varying attraction to repulsion range ratio, akin to the stripe-like patterns in 3D discussed in ref. 27, 48 and 49. Moreover, there might be a hexatic phase specific to two-dimensional systems only. An intermediate hexatic phase has been observed, e.g., in the freezing/melting of a two-dimensional dispersions of hard disks, and of disks interacting by an attractive square-well potential⁵⁰ which can be described by the celebrated Kosterlitz–Thouless–Halperin–Nelson–Young (KTHNY) theory of thermally driven dislocation or disclination pairs unbinding transitions.^51–53

Exploring the structure and phase behavior of Q2D-SALR systems is particularly relevant for colloidal (or protein) particles confined at, or close to, a planar interface. In biology, the diffusion and self-organization of phospholipids and membrane-proteins into finite-sized domains, such as “lipid rafts”^54–56 or protein clusters,^57–59 are significant features in signal transduction, membrane sorting, protein processing, and virus trafficking.^60–62 It has been evidenced that the interplay of SA (due, e.g., to lipid-mediated depletion, wetting, and hydrophobic mismatch, or direct chemical interactions between amino-acid side chains) and LR (induced, e.g., by mechanical deformation or fluctuations of the membrane) is crucial for the formation of membrane-protein clusters.^63–68 In this context, systematic investigations of the phase behavior, and the structural and cluster properties, of Q2D-SALR particles are thus indispensable. Surprisingly, despite its fundamental importance, only little effort was invested so far in exploring systematically the structure and dynamics of (Q)2D-SALR systems. There exists, however, some theoretical and simulation works on (Q)2D-SALR systems showing that particles can self-assemble into various exotic microphases either in bulk fluid,^69–72 or when trapped (pinned) by external forces^73–75 or by confinement.²⁵

In this work, using Langevin dynamics (LD) simulations, we explore the phase behavior, and the structural and clustering properties of a planar monolayer of SALR Brownian particles restricted to in-plane motion. The particles are described as Brownian spheres interacting via short-range generalized Lennard-Jones and long-range repulsive screened Coulomb forces of comparable interaction strengths. Our study covers a broad range of attraction strength and particle concentration values, and we identify and analyze the (Q)2D phases in this parameter space region. In addition to analyzing the observed cluster morphologies, we identify and quantify various (Q)2D indicators for the transition from a higher-temperature, non-clustered DF phase to a lower-temperature EC phase. The LD simulation results enable us not only to study static properties of Q2D-SALR systems as presented in this work, but also to gain (to some extent) insight into the dynamic properties of the different phases. The dynamics of Q2D-SALR systems will be discussed extensively in a subsequent paper which further addresses effects of interparticle hydrodynamic interactions not accounted for in LD simulations.⁷⁶

The paper is organized as follows. Section 2 describes the employed pair potential and the LD simulation scheme for the Q2D-SALR model system of Brownian particles. Furthermore, we recap the theoretical tools used to identify and analyze the morphology of the encountered phases. Section 3 includes the main results of the paper. Notably, a generalized phase diagram is mapped out by analyzing the shapes of the cluster size distribution function. Akin to three-dimensional SALR systems, the Q2D-SALR diagram in reduced units is subsequently amended by a metastable binodal curve for a representative short-range attractive two-dimensional square-well (SW) fluid, obtained from a second-order thermodynamic perturbation theory calculation in local compressibility approximation. Furthermore, the intermediate range order prepeak in S(q) is investigated. Various indicators for the DF to EC phase transition in Q2D systems are discussed and compared with their 3D analogues. The level of hexagonal ordering in clusters belonging to different cluster phases is carefully examined. Section 4, finally, comprises our summary with conclusions.

2 Theoretical background

We discuss here the employed Q2D-SALR model of Brownian particles and the Langevin dynamics simulation scheme, in conjunction with the theoretical tools used in analyzing the simulation data.

2.1 SALR pair potential

Our simplifying Q2D-SALR model system consists of a planar monolayer of monodisperse spherical particles embedded in a bulk fluid, with the motion of the particle centers confined to in-plane motion. The pair potential, u(r), for two particles i and j at center-to-center distance r = r_ij is thus taken to be equal to that of SALR particles in a three-dimensional dispersion, but with the constraint that r is an in-plane pair distance. We use here a standard soft SALR pair potential, namely^28,69,77The first contribution in eqn (1) is a generalized 100-50 Lennard-Jones pair potential, here referred to as u_LJ(r;ε), consisting of a hard-sphere (disk)-like steep repulsion part, and a short-range attraction part of strength ε > 0, where σ is the soft particles diameter. The second contribution, referred to as u_Y(r), is a long-range, repulsive screened Coulomb potential (Yukawa potential), where ℓ_B is the Bjerrum length of the solvent, Z_eff is the effective particle charge in units of the proton charge, and λ_D is the Debye screening length due to counter- and electrolyte ions dissolved in the solvent. The continuity of the employed SALR pair potential u(r) renders it well-suited for Langevin dynamics simulations for which the potential derivative is required.

Concerning the results shown in this paper, the non-dimensionalized Debye screening length (repulsion range) is fixed to λ_D/σ = 1.794, and the non-dimensionalized electrostatic repulsion parameter to ℓ_BZ_eff²/σ = 3.588. The non-dimensionalized short-range attraction strength, ε* = ε/k_BT, is varied in between 2–20, which amounts to a variation of the dimensionless effective temperature defined by

across the interval 0.05–0.5.

The employed SALR potential attains its minimal value, u(r_min) = u_Y(r_min) − ε ≈ 2k_BT − ε, at the pair distance r_min ≈ 1.01σ which is practically independent of the attraction strength. The attraction range parameter λ_att(ε) = r₀(ε) − σ is defined here in terms of the distance r = r₀(ε) located to the right of the potential minimum where u(r₀) = 0. This distance is given approximately by r₀(ε) ≈ (2ε*)^1/50σ for ε* ≥ 4, so that λ_att < 0.08σ in the considered range of attraction strength values. Considering that the electrostatic repulsion range, λ_D = 1.749σ, is distinctly larger than λ_att, u(r) qualifies indeed as a SALR potential. The motivation for the employed λ_D is that it is comparable to the repulsion range in previous simulation works on 3D-SALR systems.^2,23,78 Considering a colloidal suspension as an example, the employed reduced Debye screening length, λ_D/σ, is that of a dilute aqueous 1 : 1 electrolyte solution of concentration ≈3 µM, and colloidal particle diameter σ ≈ 100 nm.

2.2 (Q)2D Langevin dynamics simulations

Consider N_p disk-like and monodisperse particles interacting pairwise via the SALR pair potential given in eqn (1). The (underdamped) two-dimensional Brownian motion of the particles of mass M and diameter σ is described by the set of N_p coupled two-dimensional Langevin dynamics (LD) equations,^79,80

Mr̈_i(t) + γ₀ṙ_i(t) = −∇_iV(t) + Γ(t), (3)

for i = {1,…,N_p}. Here, r_i(t), ṙ_i(t) and r̈_i(t) are the two-dimensional position vector, velocity, and acceleration of the ith particle at time t, and γ₀ is the single-particle translational friction coefficient, associated with the single-particle translational diffusion coefficient by D₀ = k_BT/γ₀. Bearing in mind that a planar monolayer of Brownian spheres embedded in bulk fluid is considered in actuality, the three-dimensional Stokes–Einstein relation D₀ = k_BT/(3πησ) can be used for stick hydrodynamic boundary conditions, where η is the shear viscosity of the suspending fluid.

Furthermore, Γ(t) is an orthogonal, Gaussian white-noise random force vector, fully characterized by its zero mean,

〈Γ(t)〉 = 0, (4)

and its two-dimensional covariance matrix, given in dyadic notation by

〈Γ(t)Γ(t′)〉 = 2γ₀k_BTÎδ(t − t′). (5)

Here, Î is the two-dimensional unit tensor consistent with the isotropic random bombardment of the particles by the small solvent molecules. In eqn (5), the brackets 〈⋯〉 indicate a time average, k_B denotes the Boltzmann constant, and T is the temperature. The noise term acts also as a thermostat keeping the temperature of the dispersion fixed at T. The N_p Langevin equations are coupled via the direct force terms,acting on each particle i owing to the SALR pairwise interactions with the other particles. Here, V is the pairwise additive N_p-particles potential energy, r̂_ij is the unit vector pointing from the centre of particle i to that of particle j, and u′(r) is the analytically obtained derivative of the soft SALR potential given in eqn (1).

In our LD simulations, the pair potential is truncated at r = 5σ, and shifted slightly downwards to maintain its continuity everywhere.⁷⁷ Using a two-dimensional square as the primary simulation box, periodic boundary conditions and the minimum image convention are invoked. The LD equations include no hydrodynamic coupling of particle positions and orientations, i.e. solvent-mediated inter-particle hydrodynamic interactions are disregarded. Hydrodynamic coupling would turn the considered monolayer system into a genuine quasi-2D system where the lateral (out-of-plane) extension of the fluid and particles matters. We are dealing here with static properties only where the influence of hydrodynamic interactions is negligible, provided the system is in thermodynamic equilibrium.

Energy in the employed LD simulation algorithm is expressed in units of the thermal energy k_BT, length in units of σ, force in units of k_BT/σ, and time in units of which is the single-particle momentum relaxation time. The according unit for the single-particle translational friction coefficient γ₀ is .

The LD simulations integrate the Langevin equation using a small time step of approximately . The simulations run for at least 10⁸ time steps, or 2 × 10⁸ steps for attraction strengths larger than 10k_BT. This is nearly 2 × 10⁴ times the particle momentum relaxation time. Systems are prepared at different area fractions, with initial particle positions assigned according to a uniform distribution. Velocities are initialized according to a Maxwellian distribution at a temperature one order of magnitude larger than the target value. Each system is then quenched to the target temperature, and equilibration is monitored subsequently by calculating the average potential energy per particle and k_BT, as function of the LD simulation time t. We note here that methods such as simulated annealing and parallel tempering have been implemented to improve the simulation convergence towards the stationary state.^31,48,49,81 As shown by Imperio and Reatto in ref. 48 parallel tempering simulations actually achieve lower energy states, and allow a more precise computation of thermodynamic properties than plain simulations, although the overall picture on the structural features of the explored patterns does not change qualitatively.

While LD simulations allow also for the calculation of dynamic properties including time correlation functions (but with hydrodynamic interactions disregarded), in the present work only static structural properties are considered which are not influenced by inertia effects (on a coarse-grained time scale much larger than ) depending on the particle mass M. The dynamics of our Q2D-SALR model viewed on different time scales, and with hydrodynamic interaction effects included using additional multi-particle collision dynamics simulations, will be the topic of a subsequent paper.⁷⁶

2.3 Binodal of the SA reference system

Motivated by the works of Godfrin et al.²³ on the phase behavior of 3D-SALR systems, for the considered Q2D-SALR model we discuss next our calculation of the binodal of the SA reference system. We define the pair potential, u_ref(r), of the short-range attractive reference system in terms of the SALR potential u(r) in eqn (1) as

u_ref(r) = u(r)Θ(r₀(ε) − r), (7)

where Θ is the unit step function. The cut-off distance, r₀(ε), beyond which the reference potential vanishes is the second zero-value crossing point of u(r). While the attraction range λ_att(ε) = r₀(ε) − σ increases with increasing attraction strength ε, in the considered attraction strength interval, it amounts to less than 8% of the particle diameter, as noted already in Section 2.1. Thus the reference system should qualify for the Noro–Frenkel ELCS behavior, provided this extended law (verified so far for 3D-SA systems only) carries over from three to two dimensions.

To determine the gas–liquid-type binodal of the reference potential system, we employ 2nd-order (discrete) thermodynamic perturbation theory,^{22,23,82–84} but now applied to a two-dimensional fluid. To avoid cumbersome numerical calculations, triggered by the complex form of u_ref(r), we map the reference system summarily onto a two-dimensional square well (SW) fluid of comparable attraction range. The pair potential of the 2D-SW fluid is

where λ_SW − σ > 0 denotes the SW attraction range, and ε_SW > 0 the depth of the SW potential. We have identified here the hard-disk diameter of the SW fluid with the soft diameter σ of the SA reference potential, which is justifiable since the Barker–Henderson effective diameter of the reference system is only slightly different from σ (i.e., by less than 2%). For simplicity, we identify next λ_SW summarily with the cut-off distance r₀(ε) of the reference potential, averaged over the range of covered SALR attraction strength values, i.e.where ε* = ε/k_BT. The attraction range, λ_SW − σ, of the so-specified 2D-SW system amounts then to only 6% of the particle diameter. In another approximation step, we identify the SW potential depth, ε_SW, with the attraction strength ε of the SALR potential.

The binodal line of the so-specified 2D-SW fluid can be easily calculated using the 2nd-order thermodynamic perturbation theory parameterization in local compressibility approximation, as worked out by Trejos et al.⁸⁵ in the framework of the statistical associating fluid theory for potentials of variable range (SAFT-VR).^86–88 In a nutshell, thermodynamic perturbation theory computes the free energy approximately by treating the pair potential of the interacting particles system as perturbation of the hard-disk potential. For the radial distribution function (RDF) input of hard disks, Trejos et al. use an analytic superposition of the exact RDF of hard rods (Tonk gas) and the approximate RDF of a three-dimensional hard-sphere system in Percus–Yevick approximation.

The explicit form of the 2D-SW reference system binodal, as determined in this work using the semi-analytic scheme by Trejos et al., is shown in the upcoming Section 3, and it is amended by a discussion of the applicability of 2nd order perturbation theory to 2D fluids for which thermal fluctuations are stronger than in 3D.

2.4 Identification of phases

To gain insight into the phase behavior of the Q2D-SALR dispersion, using the LD simulation data we analyze the cluster size distribution (CSD) function defined aswhere N(s) is the mean fraction of particles which are members of a cluster of size s, i.e., of a cluster consisting of s particles. Here, 〈⋯〉 denotes an average over representative particle configurations, n_c(s) is the number of clusters of size s for a given configuration, and N_p is the total number of particles in the two-dimensional primary simulation box. Since sn_c(s) is the number of particles in s-clusters, it holds for each configuration that , which in turn implies the normalization condition . The mean cluster size is given by . Whether two interacting particles i and j belong to the same cluster is determined by their mutual center-to-center distance r_ij, which for a common cluster is required to be smaller than the threshold distance, r_max(ε), with r_max(ε) > r₀(ε), where the potential in eqn (1) attains its local maximum. For particles pair distances r_min < r_ij < r_max(ε), the two particles attract each other. The cluster threshold distance varies in between 1.11 < r_max(ε)/σ < 1.17 for the considered range of attraction strength values. The Q2D-SALR phases discussed in Section 3 are identified essentially on basis of the generated CSD functions.

2.5 Hexagonal order parameter

It will be shown in Section 3 that particles with SALR interactions tend in two-dimensional settings to form clusters with hexagonal ordering.⁷¹ The degree of hexagonal ordering can be analyzed, e.g., using the particle-specific local hexagonal orientational order parameter |q⁽ⁱ⁾₆|² ∈ [0,1] for a given particles configuration which quantifies the sixfold planar orientational order in the neighborhood of a particle i.^{51,52,89–92} Its complex-valued amplitude factor is given bywhere the sum runs over the six (intra- or inter-cluster) nearest neighbors of particle i. Furthermore, α_ij is the angle between the two-dimensional vector r_i − r_j of length r_ij pointing from particle j to particle i, and an arbitrarily selected in-plane axis. Note that for a cluster particle i, it is more likely that a neighboring particle j belongs to the same cluster than to another one. For |q⁽ⁱ⁾₆|² = 1, there is perfect local hexagonal ordering around a particle i, whereas |q⁽ⁱ⁾₆|² = 0 holds for (orientationally) randomly distributed particles in its vicinity.

Using our LD simulations data for the orientational order parameter subdivided into thin bins including values in between zero and one, we have furthermore calculated the orientational distribution function, P(|q₆|²), where P(|q₆|²) d|q₆|² is the probability of finding particles in the two-dimensional system with hexagonal order parameter values in an infinitesimal interval d|q₆|² centered around |q₆|². This continuous distribution function is normalized according to

For randomly distributed particles, P(|q₆|²) decays monotonically from value one at |q₆|² = 0 down to zero for |q₆|² < 1, while for significant hexagonal ordering characteristic peaks are developed in P(|q₆|²). Hexagonal ordering in the identified Q2D-SALR phases is thoroughly discussed in Section 3.3.

3 Results

We explore the microstructure and phase behavior of the employed Q2D-SALR system by considering at least 8 different values for the reduced effective temperature T* = k_BT/ε in the range [0.05–0.5], and 4 different values for the area fraction of particles, ϕ_2D = n_2Dπσ²/4, in the interval [0.1–0.5], where n_2D = N_p/L² is the area concentration of particles. The number, N_p, of particles in the LD simulations is fixed, while the areal simulation box length L is varied for setting up the considered values of ϕ_2D.

For each parameter pair (ϕ_2D,T*) refered to as state point, quantities of interest such as N(s) and P(|q₆|²) are calculated as time averages in a LD production simulation run, following an equilibration run executed over a sufficiently long time span until steady-state is reached. It was checked that the steady state value of the average energy is independent of the initial (uniform) particle configuration and the initial Maxwellian velocity distribution. However, as in ref. 48,81 a very slow logarithmic drift of the average potential energy is monitored over the entire duration of the thermalization simulation run, for cluster and percolated states at ε ≥ 10k_BT. If not specified otherwise, the presented results are for N_p = 1024 particles in the two-dimensional primary simulation box of area L × L. The particle number N_p = 1024 commonly used in our simulations is in line with previous 2D and 3D studies.^{23,29,48,49,81} A larger particle number, N_p = 4096, was used in our simulations for specific state points (ϕ_2D,T*) = (0.5, ≤0.1), for improved statistics and as a check of internal consistency.

3.1 Generalized phase diagram

As it is thoroughly discussed in the context of 3D-SALR systems,^3,23 the shape of the discrete cluster size distribution (CSD) function N(s) introduced in eqn (10) is indicative of the phase state of a considered system. For a total of 34 different Q2D-SALR systems of particles interacting by the pair potential described in Section 2.1, which are specified by respective (ϕ_2D,T*) state points, we have calculated the according N(s) functions. Recall from Section 2.4 that N(s) is the mean number of clusters consisting of s particles. Inspection of the N(s) curves has led us to the identification of four different disordered (modulated) phases described in the following.

Fig. 1(a) depicts the shapes of N(s) characteristic of the four identified (Q)2D-SALR phases. Akin to the phase diagram of 3D-SALR systems discussed in ref. 3 and 23, the identified phases encompass a dispersed fluid (DF) phase (blue squares) consisting mainly of particle monomers, an equilibrium cluster (EC) phase (red squares) with clusters having a preferential size, a random percolated (RP) gel-like phase (green diamonds) exhibiting area-spanning clusters formed basically by monomers, and a cluster percolation (CP) gel-like phase (magenta triangles) where the system-spanning clusters consist of finite cluster subunits. As noticed from the figure, the N(s) describing a DF phase system (depicted in blue) at state point (ϕ_2D,T*) = (0.2,0.5) decays monotonically with increasing s, and has its maximum at s = 1. In contrast, the CSD function of an EC fluid phase system (in red) at state point (0.2,0.125) reveals a local peak, indicative of the preferential cluster size. A RP system (depicted in green) at (0.5,0.25) is characterized by a slow, monotonic decay of N(s) terminating in a pronounced peak at s ≈ N_p, indicative of clusters spanning the two-dimensional primary simulation box with N_p = 1024 particles. Furthermore, the N(s) of a CP state system (in magenta, at (0.38,0.1)) showcases likewise a system-spanning cluster peak at s ≈ N_p, plus additional peaks reflecting smaller cluster subunits.

Fig. 1 (a) Cluster size distribution function, N(s), of Q2D-SALR systems for area fraction and effective temperature state points (ϕ_2D,T*) corresponding to the dispersed fluid (DF: blue squares), equilibrium cluster (EC: red circles), random percolation (RP: green diamonds), and cluster percolation (CP: magenta triangles) phases. (b) Corresponding reduced Q2D-SALR phase diagram. The considered 34 state points are marked according to the identified four phases. The lines are the binodals of the two-dimensional square-well (SW) reference fluid, for attraction range parameter values λ_SW = 1.075 (dashed blue), 1.06 (solid black) and 1.05 (dashed red), respectively, calculated using the 2nd order thermodynamic perturbation theory method proposed by Trejos et al.⁸⁵ All state points are normalized by the critical point values (black star) of the 2D-SW reference fluid with λ_SW/σ = 1.06. The horizontal white rectangle includes the four state points whose static structure factors are shown in Fig. 3(a).

Based on the CSD functions determined for the 34 considered state points, we obtain the Guggenheim-type (reduced) phase diagram displayed in Fig. 1(b). The (Q)2D-SALR phase diagram resembles globally that for 3D-SALR systems.²³ At comparatively high effective temperatures T* but small area fractions ϕ_2D, the systems are in the monomer-dominated DF phase space region, while the gel-like RP phase region is entered for ϕ_2D ≳ 0.37. Notice in comparison that the maximal area fraction of a planar monodisperse hard-disk layer is , attained at hexagonal closed packing. The EC and CP cluster phases are observed for lower T* ≲ 0.2 only, and for non-small ϕ_2D. The transition from the EC fluid phase to the CP gel-like phase occurs roughly for area concentrations exceeding ϕ_2D ≈ 0.35.

Akin to 3D-SALR systems, the lower-temperature two-dimensional EC and CP cluster phases in the phase diagram of Fig. 1(b) are qualitatively demarcated from the higher-temperature, non-clustered DF and RP phases by the (metastable) gas–liquid coexistence line (binodal) of a SA reference fluid of particles including short-range attractions only, with the SA pair potential defined in eqn (7) of Section 2.3. For simplicity, and as detailed in Section 2.3, we have calculated the binodal approximatively only, by mapping the SA system onto a two-dimensional square-well (SW) fluid whose attraction range parameter is fixed, according to eqn (9), to λ_SW/σ = 1.06. As another simplification, the SW potential depth ε_SW is identified with the attraction strength ε of the SALR potential. We use the 2D-SW fluid here as a simple model for the SA reference fluid, since different from the 3D case⁹³ no analytic expression is available in 2D for the free energy of the limiting case of an adhesive hard-disk system.

We have calculated the binodal characterizing the metastable gas–liquid phase coexistence of the 2D-SW fluid using second-order thermodynamic perturbation theory in local compressibility approximation, as conveniently worked out in a semi-analytic form by Trejos et al.⁸⁵ The resulting binodal for λ_SW/σ = 1.06 is shown in Fig. 1(b) as the solid black curve. Its critical point is located at . These critical point values are used to scale the area fraction ϕ_2D and effective temperature, T* = k_BT/ε, respectively, in the depicted phase diagram. On a qualitative level, the 2D-SW reference fluid binodal separates the non-clustered from the clustered phases, supporting thus the physical interpretation of the EC and CP cluster phases as modulated phases originating from gas–liquid-type demixing, macroscopically frustrated by the long-range interaction part of the SALR potential. In accord with physical expectation, and as shown in the figure, the binodal is shifted upwards (and its critical point temperature and concentration shifted to larger and smaller values, respectively) when the attraction range λ_SW is increased (see dashed lines).

We refrain from a more quantitative discussion of SA reference fluid binodals, since first these are not central to the present discussion of Q2D-SALR systems, and second since the binodals in Fig. 1(b) are mean-field results based on thermodynamic 2nd-order perturbation theory. The mean-field binodals overestimate and underestimate ϕ_c, owing to the underestimated thermal fluctuations. These fluctuations are more pronounced for smaller-range attractions as the ones considered here, and are stronger in 2D than in 3D. In this context, notice in Fig. 1(b) the strong sensitivity of the 2D mean-field binodal on the attraction range, which is likely due to the underestimation of thermal fluctuations in the employed 2nd-order perturbation scheme. Moreover, the mean-field coexistence curves underestimate the pronounced flatness of binodals generated by computer simulations,⁸⁵ and the mean-field critical point exponent β = 1/2 in the near-critical dependence of the liquid–gas concentration difference, determined along the mean-field binodal close to the critical point, differs distinctly from the exponent β = 1/8 of the 2D Ising model universallity class. The latter critical exponent is expected to apply to two-dimensional SA fluids. As noted already in Section 1, Godfrin et al.²³ discuss binodals of the 3D-SA reference fluid based on the Noro–Frenkel extended law of corresponding states (ELCS), known to hold to good accuracy for 3D fluid systems with short-range attractions where the gas–liquid type coexistence is metastable against crystallization. To our knowledge, no general validation of the ELCS has been made to date for 2D-SA fluids using computer simulations. Such a validation is more demanding in 2D than in 3D owing to stronger thermal fluctuations in 2D, and an accordingly more pronounced flatness of binodals in the (ϕ_2D,T*) diagram complicating the determination of critical points. In 2D, there is also an enhanced tendency for crystallization.

It is instructive to analyze simulation snapshots of the four detected phases. Fig. 2 includes configurational snapshots of Q2D-SALR particles for the same four state points (ϕ_2D,T*) as in Fig. 1(a). Particles belonging to clusters of different size s are distinguished by colors, selected according to the vertical color bar spanning the range from s = 1 up to s = s_max, where s_max denotes the maximal cluster size at specific (ϕ_2D,T*). In the DF system snapshot depicted in Fig. 2(a), mostly monomers are observed, but also loosely packed small clusters of varying size s are occasionally recognized. In Fig. 2(b) visualizing the EC phase of equilibrium clusters, the clusters have a preferred size (here, around s = 10), due to the increased attraction (decreased T*) in comparison with the DF phase. The RP phase snapshot in Fig. 2(c) features system-spanning clusters. It can be viewed as a crowded DF phase. Accordingly, the RP phase lacks orientational order of clusters. In contrast, the cluster percolated (CP) phase system in Fig. 2(d) has a system-spanning cluster, formed by compact and locally elongated structures with pronounced internal hexagonal order.

Fig. 2 Simulation snapshots of Q2D-SALR systems representing the four identified phases: (a) dispersed fluid (DF) phase at (ϕ_2D,T*) = (0.2,0.5), (b) equilibrium cluster (EC) phase at (ϕ_2D,T*) = (0.2,0.125), (c) random percolation (RP) phase at (ϕ_2D,T*) = (0.5,0.25) and (d) cluster percolation phase at (ϕ_2D,T*) = (0.38,0.1). The vertical cluster size color code bar ranges from s = 1 (monomers in blue) up to a maximal cluster size s = s_max at specific (ϕ_2D,T*) marked in red. The four considered area fractions correspond to simulation boxes of different sizes, since N_p = 1024 is held constant. The spatial resolution is differently selected for each of the four panels for better visibility.

A comment is in order about why globally the same phases are observed in the present (Q)2D study as in the works by Godfrin et al. on 3D-SALR systems. This can be attributed to our selection of the SALR potential parameters {ε/k_BT,l_BZ_eff²/σ,λ_D/σ} similar to those by Godfrin et al., and our focus on the phase space region around the SA fluid reference binodal. Well outside the explored phase space region, there is the possibility of finding ordered modulated phases (at low T*), as well as a hexatic phase encountered in 2D systems only. Finally, in the explored phase space region, we have no evidence of a Wigner (cluster) glass which different from the CP phase has no percolated clusters.

3.2 Prepeak characteristics at DF to EC phase transition

As mentioned in Section 1, the transition from the equilibrium DF to the modulated equilibrium EC phase of 3D-SALR systems is signalled by a characteristic height, S(q_c), of the static structure factor prepeak (intermediate range order peak) centered around the wavenumber q_c, in conjuction with a characteristic half-width, 2/ξ_T, of the prepeak where S(q_c ± 1/ξ_T) ≈ S(q_c)/2. According to Monte-Carlo simulation data by Godfrin et al.,²³ this transition occurs in 3D when the prepeak height exceeds the characteristic value S(q_c) ≈ 2.7. Moreover, as argued by Bollinger and Truskett,^33,47 for a clear-cut identification of the DF to EC transition, the prepeak height criterion should be amended by the condition that the length ξ_T slightly exceeds the Debye screening length, λ_D, of the repulsive long-range Yukawa part of the SALR potential in eqn (1). Explicitly, ξ_T ≳ λ_D is required, or more generally for an arbitrary SALR potential that 2.0 ≲ ξ_T ≲ 3.0. The length ξ_T is also refered to as thermal correlation length. Both in 3D and 2D, ξ_T is obtained from a fit of the prepeak region of S(q) around q ≈ q_c, using the inverse (Bollinger–Truskett) quadratic formWhile this form bears some similarity with the Ornstein–Zernike expression of the structure factor near zero wavenumber for a near-critical fluid, the Bollinger–Truskett form is obviously not even in q for q_c > 0.

Fig. 3(b) shows LD simulation results for the static structure factor,

of our Q2D-SALR model at fixed temperature T* = 0.167 (i.e., ), and for four different area fractions ϕ_2D = 0.1–0.5, representative of three out of the four detected phases. The brackets indicate an average over a stationary ensemble of particle configurations, and q is the scattering wavevector of magnitude q ≥ 2π/L (with L the linear dimension of the simulation square), selected according to the invoked periodic boundary conditions. The prepeaks visible at wavenumber position q_cσ ∼ 1.5 increase in height and become more narrow with increasing area fraction. A distinctive prepeak (in blue) is visible already for the DF phase system at ϕ_2D = 0.1. While q_c remains nearly constant for the DF and EC systems, it is shifted to a somewhat smaller value for the CP system (magenta curve) at ϕ_2D = 0.5, for which the nearest-neighbor shell peak of S(q) at q_m ≈ 2π/σ > q_c is distinctly higher than the prepeak. That q_c is rather insensitive to changes in particle concentrations is observed also for 3D-SALR systems. Akin to the 3D case, the concentration insensitivity of q_c can be qualitatively explained using a simple random phase approximation where the static structure factor of the (Q)2D-SALR system is related approximately to the static structure factor, S_hd(q), of a hard-disk reference system. On noting that S_hd(q ≲ q_c) ≈ S_hd(q = 0), it is predicted in random phase approximation that q_c is independent of ϕ_2D, and determined only by the competing attraction/repulsion interaction parameters of the SALR potential. In actuality, q_c can be mildly ϕ_2D-dependent for fixed pair potential.

Fig. 3 (a) Static structure factor, S(q), as function of the non-dimensionalized wavenumber qσ, for reduced temperature T* = 0.167 (ε = 6k_BT) and four different area fractions as indicated. The according four state points are included in the horizontal white rectangle in Fig. 1(b). The dashed vertical arrow points to the locations, q_cσ, of the intermediate range order prepeaks. (b) Static structure factor prepeak height S(q_c) (logarithmic left ordinate, filled symbols) and reduced thermal correlation length ξ_T/σ (right ordinate, open symbols) as functions of area fraction. The color code for the different phases is as in Fig. 1. The solid (dashed) horizontal line marks the identified transition value S(q_c) ≈ 1.4 (ξ_T/σ ≈ 2), separating approximately the DF from the EC phase for ϕ_2D ≲ 0.4.

Fig. 3(b) summarizes our results for the height, S(q_c), of the prepeak of S(q) (filled symbols and left ordinate), and for the corresponding thermal correlation length ξ_T/σ (open symbols and right ordinate), the latter obtained from a best fit of the simulated (Q)2d-SALR prepeak using the analytic form in eqn (13). Different from the figure part (a), T* is not fixed in figure part (b) so that all explored state points are considered, the colors labeling different phases are selected as in Fig. 1.

Quite interestingly, the characteristic prepeak height S(q_c) ≈ 1.4 signaling the transition from the DF to the EC phase for ϕ_2D < 0.4 (solid horizontal line) is distinctly smaller than the according characteristic value S(q_c) ≈ 2.7 in the 3D. This is in accord with our expectation that cluster formation is facilitated in 2D, since the local order in 2D is sixfold in all observed phases. As noted in Section 1 already, the DF to EC transition height S(q_c) ≈ 2.7 in 3D-SALR systems²³ is close to the Hansen–Verlet freezing transition criterion values, S(q_m) ≈ 2.85–3.1, for the nearest-neighbor peak at q_m. The latter peak heights signal a macroscopic fluid to crystal (freezing) transition in 3D systems with dominantly repulsive pair interactions. The freezing transition value S(q_m) ≈ 2.85 is attained for hard spheres at volume fraction ϕ_3D = 0.494, while the upper value 3.1 is characteristic for the freezing of complex fluids having long-range repulsive interactions such as charge-stabilized suspensions of low salt ion content. This should be contrasted with the according Hansen–Verlet freezing transition peak values for 2D systems, obtained from simulations, which span the range from S(q_m) ≈ 4.4 for a hard-disk system at freezing for ϕ_2D ≈ 0.70, up to S(q_m) ≈ 5.7 for particles repelling each other by a soft 1/r¹² pair potential where freezing occurs at ϕ_2D ≈ 0.77.^94–97 For super-paramagnetic particles at a planar gas–liquid interface interacting by a long-range repulsive 1/r³ pair potential, an even larger Hansen–Verlet value S(q_m) ≈ 10 is found experimentally at the hexatic to crystalline phase transition point.⁹⁸ Characteristic of 2D systems at freezing is a strong change of S(q_m) across the freezing point.^97,98

The freezing transition peak height, S(q_m), of 2D systems with (dominantly) repulsive interactions is thus more than three times larger than the prepeak value S(q_c) ≈ 1.4 signalling DF to EC transition in (Q)2D-SALR systems. However, one needs to distinguish the different mechanisms driving the respective transitions. Macroscopic freezing is driven by particles packing/caging effects. While the mesoscale DF to EC transition is driven by the competition of SA and LR pair forces. In this context, notice from Fig. 3(b) the pronounced increase in S(q_c) with increasing area (packing) fraction, causing for fixed temperature the subsequent transitions DF → EC → CP. The distinctly larger Hansen–Verlet peak heights S(q_m) as compared to 3D systems, and the larger area concentration values at 2D freezing, can be attributed to stronger critical thermal fluctuations in 2D, which necessitates stronger interparticle correlations to stabilize two-dimensional crystals exhibiting quasi long-range order only.

Having identified the prepeak height S(q_c) ≈ 1.4 as an approximative indicator of the DF to EC transition in (Q)2D-SALR systems, we analyze next our simulation results for the thermal correlation length ξ_T depicted in Fig. 3(b) (open symbols and right ordinate). As indicated by the horizontal dashed line, the EC phase region for ϕ_2D ≲ 0.4 is reached approximately, in progressing from the DF phase region, provided ξ_T ≳ 2σ. This threshold value implies that to form (Q)2D equilibrium clusters, ξ_T should exceed the Debye screening length λ_D in the LR Yukawa part of the SALR potential in eqn (1), whose range is fixed to λ_D = 1.794σ in our study. Our finding for the characteristic value of ξ_T is therefore similar to the one by Bollinger and Truskett for 3D-SALR systems,³³ but we cautionally recall that no span of different λ_D values is covered in the present simulations.

It is instructive to analyze asymptotically the two-dimensional, real-space Fourier back transform,

of the Bollinger–Truskett expression in eqn (13) describing the prepeak region of S(q) around q_c. Here, J₀(x) is the zeroth-order Bessel function of first kind. From a large-distance r analysis of the integral based on the residue calculus, we obtainwhere B = (q_c² + 1/ξ_T²)^1/2 and ψ = arctan(1/q_cξ_T). The oscillatory first term in the asymptotic expression for the real-space prepeak total correlation function h_c(r), valid for r ≫ ξ_T, expresses undulations in the particle pair correlations on length scale 1/q_c, while its envelope decays in two dimensions as . This clarifies the physical meaning of ξ_T as the range of pair correlations. The second polynomial term decaying proportional to 1/r³ arises from the linear in q contribution of the denominator, 1 + ξ_T²(q − q_c)², of the integrand, and it dominates the asymptotics of h_c(r) at large distances. The second polynomial term on the right-hand side of eqn (16) is non-zero, and the first term is oscillating away from a Fisher–Widom line where q_c = 0.⁹⁹ In comparison, the long-distance oscillatory part of h_c(r) in 3D decays more strongly as exp{−r/ξ_T}/r, and the polynomial term dominating the long-distance decay of h_c(r) in 3D is proportional to 1/r⁴. Our asymptotic 2D pair correlation expression in eqn (16) which is distinctly different from the 3D case, could prove useful for analyzing experimentally determined or simulated pair correlations in (Q)2D cluster phases such as studied, e.g., in ref. 32.

In concluding this subsection, akin to corresponding 3D criteria by Godfrin et al.,²³ and Bollinger et al.,^33,47 we propose the following characteristic prepeak values to localize approximately the DF to EC transition in Q2D-SALR systems: first, the prepeak height should surpass S(q_c) ≈ 1.4, and second, the reduced thermal correlation length ξ_T/σ should be not smaller than 2.0. We caution that the characteristic values for S(q_c) and ξ_T/σ indicating the DF to EC transition in 2D, may depend to some extent on the shape of the SALR potential. In 3D, the corresponding characteristic values appear to be insensitive to the SALR potential shape, which can be attributed to some extent to the validity of the Noro–Frenkel extended law of corresponding states for 3D particles systems with short-range attraction (SA). Whether this extended law remains valid in 2D – SA systems is the topic of an ongoing simulation study, into which the present authors are involved.

3.3 Local hexagonal ordering and cluster properties

In the following, we discuss first the local hexagonal ordering of (Q)2D-SALR particles, which in two dimensions is energetically favored, and progress next to the analysis of cluster properties in the different phases.

3.3.1 Radial distribution function. The local hexagonal order visible for the EC and CP cluster phases in the simulation snapshots of Fig. 2, can be quantified, e.g., by the radial distribution function (RDF) g(r), related to the static structure factor S(q) in eqn (14) by the inverse two-dimensional Fourier transformThe RDF g(r) is the conditional probability of finding a particle a center-to-center distance r apart from a given one, irrespective of whether the two particles belong to the same or different clusters, or to no cluster at all. In the LD simulations, g(r) is determined for r < L from a discretization of the average,over a stationary ensemble of statistically independent particle configurations, where L is the simulation box length, r = |r| and δ⁽²⁾(r) denotes the two-dimensional delta function. Self terms l = j are excluded from the double sum.

Fig. 4 displays LD simulation results for g(r) of clustered Q2D-SALR systems in the EC phase (upper panel, in red) and in the CP phase (lower panel, in magenta). Two state points (ϕ_2D,T*) are considered for each phase, of values given in the figure. The (dashed) curves of g(r) for the lower-concentrated system in each phase are shifted upwards by two units, to render them better distinguishable from the g(r) of the more concentrated systems (solid curves). The vertical dotted lines at mark the next, overnext and so on interparticle distances, of ascending magnitude, of close-packed disks of diameter σ centered at the vertices of an ideal two-dimensional hexagonal (triangular) lattice. These interparticle distances are readily obtained as the vector lengths of linear superpositions of the hexagonal lattice basis vectors b₁ = x̂ and , with inner product b₁·b₂ = cos(60°), where x̂ and ŷ are unit vectors along the orthogonal x and y directions, respectively.

Fig. 4 Radial distribution function, g(r), as function of reduced particle pair distance r/σ. The state points (ϕ_2D,T*) considered in the upper panel are for two equilibrium cluster (EC), and in the lower panel for two cluster percolated (CP) systems, respectively. The dashed curves are shifted upwards by two units for better visibility. Vertical dotted lines mark the locations of real-space inter-vertices distances of an ideally hexagonal 2D crystal. See main text for details.

Consider first the two EC systems in the upper panel at equal effective temperature T* = 0.063. The EC g(r) for the smaller area fraction ϕ_2D = 0.1 (dashed red curve, shifted upwards by two units) reveals already well-developed, longer-range inter-particle correlations, quantified by several sharp peaks in g(r) of sizeable heights, persisting to inter-particle distances r of about five times the particle diameter σ. The peaks in g(r) are centered at the above noted interparticle distances separating the vertices of a 2D hexagonal lattice, which indicates pronounced local hexagonal ordering inside the equilibrium clusters. These characteristic features of g(r) become sharpened and more pronounced for the EC system at the larger area fraction ϕ_2D = 0.2 (solid red curve).

Notice, in particular, the existence of a broad and very shallow peak region of g(r) with maximum at about r ≈ 6σ (pointed at by the short black arrow in the upper panel of Fig. 4), which is close to the distance r_c = 2π/q_c associated with the mean distance of neighboring clusters. Incidentally, such an extended shallow peak region of g(r) at larger pair distances is also observed in multiparticle collision dynamics simulation results for 3D-SALR systems in the EC phase region.⁴²

Consider next the two RDFs of the percolated CP systems in the lower panel of Fig. 4. These are at a temperature T* = 0.063 higher than the one of the EC systems in the upper panel, but the two area concentrations are now distinctly larger. As seen in the figure, the RDFs of the cluster percolated particles are similar to those of the EC systems, with peaks likewise positioned at the inter-vertices distances of an ideally hexagonal lattice. However, the peak heights for r/σ > 1 are in general smaller. Moreover, the shallow peak region at larger r is flattened out in comparison to the non-percolated EC systems, indicative of a more random distribution of particles at larger pair distances r ≳ 6σ (see black arrow).

3.3.2 Hexagonal ordering of clusters. To measure the degree of local hexagonal ordering of Q2D-SALR particles in the different phases, we analyze first the particle-specific local hexagonal order parameter |q⁽ⁱ⁾₆|², with its complex-valued amplitude q⁽ⁱ⁾₆ defined in eqn (11), in its dependence on concentration and effective temperature. As explained in Section 2.5, this parameter measures the local hexagonally ordered nearest neighbor particles around a particle i. It attains values in between |q⁽ⁱ⁾₆|² = 0 for orientationally disordered neighboring particles, and |q⁽ⁱ⁾₆|² = 1 for perfect hexagonal ordering of neighboring particles.

Typical configurational simulation snapshots for 20 different state points (ϕ_2D,T*) encompassing clustered EC and CP systems and non-clustered DF and RP systems, are displayed in Fig. 5. The local hexagonal order parameter of each individual particle is color-coded by the vertical bar on the right-hand side of the figure. This color coding of different hexagonal ordering should be distinguished from the four-colors code used throughout this work to label the four different phases. In the dispersed fluid (DF) snapshot for ϕ_2D = 0.2 and T* = 0.167, enclosed in a thick blue frame, for most particles it is |q⁽ⁱ⁾₆|² ≈ 0, i.e., there is no local hexagonal order. This is the case also for the random percolated (RP) system at ϕ_2D = 0.5 and T* = 0.25, whose snapshot is depicted in the top left corner of the figure enclosed by a green frame. It has no pronounced orientational order even though it has system-spanning clusters. In contrast, there is significant local hexagonal order in the 11 snapshots of the equilibrium cluster (EC) systems, shown in the figure for three area fractions and six different effective temperature values. The EC systems are the ones inside the red frames. In comparison with the DF system, there is a noticeable reduction in the fraction of particles with zero local order parameter. For the lowest-concentrated EC system at (0.1,0.125), only few hexagonally packed small clusters are formed, as reflected by a small number of particles with |q⁽ⁱ⁾₆|² ≈ 1. In fact, there is still a substantial fraction of non-clustered particles characterized by |q⁽ⁱ⁾₆|² ≈ 0. Intermediate values 0 < |q⁽ⁱ⁾₆|² < 1 characterize cluster particles with less than 6 nearest neighbors. Recall from Section 2.5 that for a cluster particle, it is more likely that its neighbors belong to its own cluster than to other ones.

Fig. 5 Typical configurational snapshots at four different area fractions, ϕ_2D, and six different effective temperatures, T*, as indicated, except for the snapshot at the top left corner, demarcated by thick green lines, which corresponds to (ϕ_2D,T*) = (0.5,0.25). The particle-specific hexagonal order parameter |q⁽ⁱ⁾₆|² ∈ [0,1] is color-coded by the vertical bar on the right hand side. Snapshots of systems representing the same phase are encompassed by a frame with color as used in Fig. 1, i.e., DF (blue), EC (red), RP (green), and CP (magenta), respectively. A selection of four double-strands hexagonal clusters are marked by small red rectangles. See main text for details.

In addressing EC systems with stronger short-range attraction where T* ≲ 0.1, an interesting change in the cluster morphology is observed. Decreasing T* from 0.125 down to 0.1, the cluster shapes progressively change from mainly disk-like to elongated, whereas the cluster sizes s have not changed significantly yet. Only for T* > 0.1 and the larger area fractions ϕ_2D = 0.1 and 0.2, a larger fraction of elongated clusters emerges. These clusters are mainly double-stranded, i.e., they consist of two strings of particles in a locally hexagonal arrangement. Four examples of these double-stranded hexagonal clusters are marked in Fig. 5 by red rectangles.

Double-stranded hexagonal particles arrangements are omnipresent in the cluster percolated (CP) phase. According snapshots of eleven CP systems are enclosed in Fig. 5 by a magenta-colored frame. As attraction is increased by lowering T*, the initially quite open percolated cluster network increasingly develops branch points, which are more numerous for the largest considered area fraction ϕ_2D = 0.5. For the lowest temperature T* = 0.167, the system-spanning network consists to a larger degree of particles belonging to perfectly hexagonally ordered cluster segments where |q⁽ⁱ⁾₆|² ≈ 1. Particles in the double-stranded cluster segments of the percolated cluster network are typically bonded to 4 neighbors. The mean length of the segments decreases with decreasing temperature.

An analogous formation of elongated clusters has been observed previously for 3D systems under two specific conditions. First, Toledano et al.⁸¹ found elongated clusters at large particle concentrations close to and above the percolation threshold, for a long-range screening length of λ_D = 2σ. Second, for distinctly shorter screening lengths of a fraction of σ, quasi-one-dimensional clustering and the formation of peculiar Bernal spirals was observed both in simulations¹⁰⁰ and experiments.⁵ The double-stranded hexagonal clusters seen in the Q2D-SALR systems of this work are in line with according findings in ref. 81. Their existence can be ascribed to inter-cluster interactions favoring anisotropic structures instead of isotropic ones to lower potential energy. We point out that enlarged short-range attraction induces anisotropic cluster formation even for small particle concentrations, far below the percolation threshold value. We argue that the outer cluster particles preferentially arrange themselves into elongated shapes, to counterbalance the net repulsion by the inner cluster particles, thus minimizing the system free energy. Especially at lower T*, it is energetically more favorable for additional particles to connect to a cluster along the elongated direction, since for a connection along the short cluster axis a larger kinetic barrier has to be surpassed.^77,101

Additional information on the hexagonal ordering of clusters in different phases is contained in the normalized, continuous probability distribution function, P(|q₆|²), of the hexagonal order parameter |q₆|² introduced in Section 2.5. Simulation results for P(|q₆|²) are presented in Fig. 6, for area fractions ϕ_2D = 0.1, 0.2, 0.38 and 0.5, respectively, and three effective temperatures as indicated. The hexagonal order distribution function of DF systems at high temperatures (cf., the blue curves in Fig. 6(a)–(c)) decays monotonically, expressesing lack of hexagonal ordering.

Fig. 6 Probability density function of the hexagonal order parameter, P(|q₆|²), for area fractions (a) ϕ_2D = 0.1, (b) ϕ_2D = 0.2, (c) ϕ_2D = 0.38, and (d) ϕ_2D = 0.5, respectively, and for three reduced temperature values T* as indicated. The color code indicating different phases is as in Fig. 1. For better visibility, some curves are dashed or semi-transparent.

A monotonic decay of P(|q₆|²) is likewise observed for the RP system (green curve) in Fig. 6(d), at state point (ϕ_2D,T*) = (0.5,0.25), since different from EC and CP systems there are basically no double-stranded hexagonal cluster (segments) in RP phase systems (see Fig. 2(c)). Akin to DF systems, the particles in the higher-temperature RP phase region forming system-spanning networks are locally mobile which allows them to readjust their relative orientations.

In contrast, at lower temperatures, the presence of hexagonally ordered clusters in the EC phase systems at ϕ_2D = 0.1 and 0.2 (red curves in Fig. 6(a) and (b)) is reflected by multiple sharp peaks in P(|q₆|²). The smaller peak at |q₆|² = 1 is due to a modest amount of six-bonds neighboring particles inside disk-shaped clusters, as it is noticed from Fig. 5. In addition, more pronounced peaks of P(|q₆|²) are located at |q₆|² = 1/9, 1/4, and 4/9. These peaks are mainly due to particles positioned at the fringe of the clusters which have less than six next neighbors. The peak at |q₆|² = 1/9, e.g., is mainly due to particles at the lateral corners of double-stranded clusters with two nearest neighbors only. Increasing the attraction strength or area fraction enhances the peak heights (dashed orange curves), since more double-stranded EC clusters of elongated shape are formed wherein the cluster particles have four or two nearest neighbors. There is an according suppression of the peak at |q₆|² = 1. Interestingly, for equilibrium cluster (EC) systems, there is a certain amount of s = 5 and s = 6 clusters with hexagonal order such that |q₆|² = 0.233 and |q₆|² = 0.198, respectively. Since these two values are close to each other, they both contribute to essentially a single peak in P(|q₆²|), visible in the figure at |q₆|² ≈ 0.2. In the percolated RP and CP phases, the clusters are very extended and hence this additional peak is not present any more.

Consider next P(|q₆|²) of the cluster percolated (CP) systems (solid and dashed magenta curves) in Fig. 6(c) and (d), where the system-spanning, branched network is put together mainly by double-stranded hexagonal cluster segments. As seen, the peak at |q₆|² = 1/9 becomes strongly enhanced with increasing attraction, consistent with an overall enthalpic strengthening of local hexagonal order for decreasing effective temperature. An enhancement of this peak with increasing concentration at fixed temperature (such as for the EC systems), could be expected on first sight also for the CP systems in Fig. 6(c) and (d). Conversely, however, this peak decreases when ϕ_2D is increased from 0.38 to 0.5. We attribute this to the increased branching of the system-spanning network, as it is visible in Fig. 2, with an accompanying loss of network particles having two nearest neighbors.

Fig. 7 shows the configuration-averaged, hexagonal order parameter, 〈|q₆|²〉, i.e., the first moment of P(|q₆|²), plotted as function of T*, for four different area fractions as indicated. The phase assignment of state points is marked by symbols and colors selected as in Fig. 1. State points of equal ϕ_2D are connected by black line segments. Notice that particles belonging to different clusters are not clearly discriminated in computing 〈|q₆|²〉, since it invokes an ensemble average over all particles. The depicted average orientational order parameter increases with decreasing temperature (T* < 0.18) towards a plateau value, whereas for ϕ_2D = 0.38 and in particular for ϕ_2D = 0.5 a sharp non-monotonic variation of 〈|q₆|²〉 is observed in the interval T* = 0.1–0.2. In comparing the cluster percolated (CP) states at (ϕ_2D = 0.5,T* ≤ 0.167) and (ϕ_2D = 0.38,T* ≤ 0.125), the latter ones have more particles in double-stranded cluster segments connecting the branch points of the system-spanning clusters. This implies less overall hexagonal ordering for the lower-concentrated CP systems, as seen in Fig. 5.

Fig. 7 Configuration-averaged hexagonal order parameter, 〈|q₆|²〉, as function of the effective temperature T*, for four different area fractions as indicated. The phase state membership of state points is indicated by symbols and colors as in Fig. 1. State points of equal ϕ_2D are connected by black line segments.

3.3.3 Mean cluster size 〈s〉. An integral measure characterizing clusters in different phases is the mean cluster size 〈s〉, i.e., the first moment of the CSD function N(s) defined in eqn (10) of Section 2.4, and displayed for different phases in Fig. 1(a).

The temperature dependence of 〈s〉 is analyzed in Fig. 8(a) for the area fractions ϕ_2D = 0.1, 0.2, 0.38 and 0.5. As expected, 〈s〉 decreases with increasing T* at constant ϕ_2D, and it increases with increasing ϕ_2D for constant T*, which is characteristic for the percolated cluster phase. In comparison, while still being large, the mean cluster size of the random percolated systems (green diamonds) at ϕ_2D = 0.5 is smaller, which conforms with the overall flat shape of the RP phase CSD function in Fig. 1(a).

Fig. 8 (a) Mean cluster size, 〈s〉, and (b) coordination number (average number of bonds per particle), 〈z_b〉, as functions of effective temperature T* = k_BT/ε. Simulation points of equal ϕ_2D are linked by black line segments. The arrows point towards increasing area fraction values ϕ_2D = 0.1, 0.2, 0.38, and 0.5. (c) Probability distribution of bond number per particle, P(z_b), for cluster percolated (CP) phase systems at five different indicated state points (magenta symbols, linked by dashed lines). (d) Coordination number, 〈z_b〉, as function of average cluster size 〈s〉, for various (ϕ_2D,T*) state points. The dashed horizontal black line marks the threshold value 〈z_b〉 ≈ 1.6 separating DF from EC phase systems. The thick solid black curve represents the empirical relation in eqn (19), as discussed in the main text. Symbols and colors indicating different phases are selected as in Fig. 1.

Interestingly, DF systems (blue squares) at higher concentrations can have larger mean cluster sizes than EC phase systems (red circles) at lower concentration and temperature, owing to packing effects coming into play at higher particle concentration. An example in case is the DF system at (ϕ_2D,T*) = (0.38,0.25) if compared with the EC system at (0.2,0.125). However, clusters of the EC phase are distinctly more persistent than the transiently formed small clusters in the DF phase. The depicted curves of 〈s〉 as functions of T* have sigmoidal shapes, reaching a plateau at low T* whose height is larger for larger ϕ_2D. A close inspection of N(s) for T* < 0.1 (data not shown here) reveals that the CSD function becomes insensitive to T* in the low temperature range, while clusters of size close to the mean value 〈s〉 are more likely to occur. In this context, note that in their 3D-SALR systems analysis based on classical nucleation theory, Zhang et al.¹⁰¹ predict an initial increase of 〈s〉 when the attraction strength is increased, in accord with our simulation results for ε/k_BT ≤ 10, while for stronger attraction (ε ≈ 16k_BT) they detect a decrease in 〈s〉.

3.3.4 Coordination number 〈z_b〉. Additional insight into the local order of clusters is gained from exploring the mean value, 〈z_b〉, of the bond number z_b, commonly refered to as coordination number. For a given configuration, z_b is the number of bonds a particle has with its neighboring particles. In two dimensions, z_b ≤ 6 owing to the SA pair potential part. As discussed in Section 2.4 in relation to N(s), a bond between two particles is formed when their center-to-center distance r is less than the distance r_max(ε) where the SALR potential has its local maximum. Experimentally, 〈z_b〉 can be determined using optical microscopy.

Simulation results for 〈z_b〉 are shown in Fig. 8(b) as functions of T*, for area fractions as in Fig. 8(a). The coordination number curves follow the same trends as the mean cluster size curves in Fig. 8(a). In particular, they share similar sigmoidal shapes, with 〈z_b〉 attaining a constant value for low temperatures. Most notably, the coordination number of the CP systems (magenta circles) at low T* is close to 4.

To elucidate this observation, we analyze next the full statistical information on the bond number z_b encoded in the discrete distribution function, P(z_b), of bond numbers, which is normalized such that . Simulation results of P(z_b) are depicted in Fig. 8(c) for five CP systems (magenta color) whose respective state points are given in the figure. For each CP system (state point), the values of P(z_b) at the discrete bond numbers z_b = 0–6 are linked by straight line segments, for better visibility. The bond number distribution functions of all five CP systems have the maximum located at z_b = 4, and the value of the maximum increases with decreasing temperature for constant ϕ_2D. The maximum of P(z_b) at z_b = 4 is due to the composition of the system-spanning network by double-stranded cluster segments, whose constituting particles are surrounded mainly by four next neighbors. While the probability for encounting monomers (i.e., for z_b = 0) is close to zero, all considered CP systems have a small but non-zero probability for z_b = 6. The CP system at state point (ϕ_2D,T*) = (0.5,0.167) falls out of line in having a 20% probability for 〈z_b〉 = 6. This enlarged probability is compensated by an accordingly smaller probability for z_b = 4 such that the normalization of P(z_b) is conserved.

Analogous to earlier findings for 3D-SALR systems,^23,100 and as suggested by the similarity of the sigmoidal curves for 〈s〉(T*) and 〈z_b〉(T*) in Fig. 8(a) and (b), respectively, a monotonic relation between 〈z_b〉 and 〈s〉 is expected also for Q2D-SALR systems. A qualitative reasoning for such a relation to exist has been given for the 3D case by Godfrin et al.:²³ with increasing mean cluster size 〈s〉, more particles are integrated into clusters which enlarges the mean number of nearest neighbors. Moreover, and once again similar to the 3D case, for a given 〈s〉 the coordination number of clustered phases can be expected to be larger than that of non-clustered ones. That this reasoning applies also to Q2D-SALR systems is evident from Fig. 8(d), encompassing all state points considered in this work. Additionally shown in the figure as solid black line is the curve of the empirical relation²³

〈z_b〉 = 1.5(ln〈s〉)^1/2, (19)

which is seen to separate the more coordinated clustered EC (red) and CP (magenta) systems from the less coordinated non-clustered DF (blue) and RP (green) systems. Interestingly enough, the very same empirical relation was shown by Godfrin et al.²³ to separate clustered systems from unclustered ones in the 3D case. In Q2D, however, the scattering of state points around this relation is distinctly less pronounced. The DF and RP points, in particular, follow the empirical relation quite closely on the depicted semi-logarithmic scale. The horizontal dashed line in Fig. 8(d) marks the mean coordination number threshold value 〈z_b〉 ≈ 1.6, separating in Q2D the DF from the clustered EC systems (in 3D, this threshold value is about 2.4.²³). The mean coordination number 〈z_b〉 ≈ 1.6 in Q2D corresponds to a critical mean cluster size of about 〈s〉 ≈ 3.1 which needs to be exceeded before equilibrium cluster formation sets in. In 3D, the critical mean cluster size based on 〈z_b〉 ≈ 2.4 and eqn (19) is 〈s〉 ≈ 12.9.

The discussed DF to EC transition indicator values in 2D and 3D are summarized and compared in Table 1. Incidentally, as shown for 3D-SALR systems,²⁸ there exists also a dynamic indicator for the DF to EC transition in terms of the mean cluster lifetime, whose slope as function of the attraction strength ε changes abruptly at the transition value. Whether this dynamic indicator remains valid also for Q2D-SALR systems will be analyzed in subsequent work.⁷⁶

Table 1 Comparison of characteristic indicator values for the DF to EC transition of SALR systems in 2D and 3D, respectively. Notice further the more general DF to cluster phases transitions indicator in eqn (19), valid both in 2D and 3D

	S(qc)	ξ T/σ	〈zb〉c	〈s〉c
2D	1.4	>2.0	1.6	3.1
3D	2.7	>2.0 ∧ <3.0	2.4	12.9

---

In a nutshell, we have established a comprehensive and coherent picture of cluster properties in Q2D-SALR fluids of pair potential form given in eqn (1), for a broad range of effective temperatures and concentrations. Using LD simulations, we have shown that Q2D-SALR clusters have pronounced local hexagonal order. For sufficiently strong short-range attraction, the system tends to form elongated double-stranded hexagonal clusters even at low area fractions. At higher concentrations, branched percolated networks are formed whose morphology is still dominated by the double-stranded hexagonal pattern of cluster segments. We caution that SALR potentials different from eqn (1) may lead to distinctly different (Q)2D mean cluster sizes and coordination numbers, for comparable attraction and repulsion strength parameters.

4 Conclusions

Using LD simulations, we have investigated the phase behavior and structural properties of Q2D-SALR systems of monolayers of Brownian spheres, confined to translational in-plane motion in an embedding fluid, and interacting by a generalized Lennard-Jones–Yukawa pair potential with competing short-range attractive (SA) and long-range repulsive (LR) interaction parts. A broad range of area concentrations and effective temperatures has been covered in our study. A generalized phase diagram was mapped out, by analysing the cluster size distribution function. We have identified four distinct phases: a dispersed fluid (DF), equilibrium cluster (EC), random percolation (RP), and cluster percolation (CP) phase. While these phases resemble those for 3D-SALR systems, they are distinctly different in their microstructure invoking different levels of hexagonal ordering.

Using a second-order thermodynamic perturbation scheme applied to a two-dimensional reference SW fluid related to the SA part of the SALR potential, we determined a binodal curve delineating approximately the non-clustered from the cluster phases. The strong sensitivity of the binodal on the attraction range, noticed from Fig. 1(b), is likely an artifact of the perturbation scheme which in 2D is distinctly less accurate close to a critical point than in 3D. For this reason, we have refrained from studying the extended law of corresponding states (ELCS) behavior of the 2D reference fluid using perturbation theory. However, a detailed 2D Monte-Carlo simulation study of the ELCS behavior of short-range attractive SW and Mie potential fluids is in progress, in which part of the present authors are involved.

Furthermore, various indicators have been analyzed for the transition from the higher (effective) temperature DF to the lower (effective) temperature EC phase region, in comparison with their 3D analogues. The DF to EC transition is signalled by the threshold value S(q_c) ≈ 1.4 of the low-wavenumber peak of the static structure factor, and the thermal correlation length ξ_T ≳ 2σ exceeding in particular the Debye screening length λ_D quantifying the repulsion range of the SALR potential. While the (Q)2D criterion for ξ_T agrees practically with that for 3D-SALR systems, the characteristic prepeak height S(q_c) ≈ 1.4 is distinctly smaller than that in three dimensions (S(q_c) ≈ 2.7). The smaller prepeak height complies with the favored local hexagonal ordering in the 2D phases, which facilitates cluster formation. Different from the 3D case, however, the Hansen–Verlet concentration peak height S(q_m) ≈ 4.4–5.7 of the static structure factor signalling freezing into a 2D crystalline phase is more than three times larger than the prepeak height at SALR equilibrium cluster formation. This can be attributed to stronger thermal concentration fluctuations in 2D destabilizing crystal formation.

We have conducted a thorough analysis of local hexagonal ordering in the four detected phases. As the attraction strength ε is increased (typically up to more than 10k_BT) at lower area concentrations, the finite-sized equilibrium clusters present in the EC phase change from nearly disk-like to elongated shapes. The observed elongated clusters are mainly double strands, i.e., they are composed of two parallel linear strings of particles arranged hexagonally. This type of 2D clusters resembles the Bernal spirals observed in 3D systems.¹⁰⁰ Bernal spiral formation is either triggered by a short electrostatic screening length, or by intra- and inter-cluster interactions in the case of larger screening lengths. The double-strands pattern is dominant also in the CP phase region of large concentrations, in form of the cluster segments of the branched percolation network.

The mean cluster size, 〈s〉, increases with increasing concentration, ϕ_2D, and increasing attraction strength (decreasing effective temperature T* = k_BT/ε), as expected. It has a sigmoidal-shape as function of increasing attraction strength, i.e., 〈s〉 grows gradually for weak attractions, followed by a sharp increase for intermediate attractions, and finally reaching a saturation value for highly attractive particles. The coordination number, 〈z_b〉, follows a similar trend as 〈s〉 for varying attraction strength. For the Q2D-SALR systems, we observe a threshold value 〈z_b〉 ≈ 1.6 separating the DF from the EC phase region. The corresponding threshold value in 3D-SALR systems is 〈z_b〉 = 2.4, which complies with the rigidity percolation threshold value at 3D covalent glass formation.^23,102

For the EC and CP cluster phases, the probability distribution of bond numbers, P(z_b), attains its maximum in the bond number range z_b ≥ 2. Provided the cluster morphology is dominated by double-stranded clusters, at concentrations below or above the percolation threshold, the peak in P(z_b) persistently occurs at z_b = 4. Intriguingly, the empirical relation by Godfrin et al.²³ between 〈z_b〉 and 〈s〉 in eqn (19), which delineates cluster from non-clustered SALR phases, remains applicable also to Q2D systems. According to this relation, the DF to EC transition threshold value 〈z_b〉 = 1.6 corresponds to a mean cluster size of 〈s〉 ≈ 3.1.

The presented simulation study of the phase behavior and structure of a generic Q2D-model of SALR particles allows to gain qualitative insight into the clustering properties of more complex Q2D systems such as particles confined to a fluid–fluid interface with viscosity contrast,¹⁰³ and into biologically relevant membrane protein aggregation effects at the intracellular level. A natural extension of this study is the exploration of the (cluster) dynamics of the employed Q2D-SALR model, including the role of the hydrodynamic interactions for diffusion and surface-rheological properties. Work in this direction is in progress, and its results will be reported in a subsequent paper.⁷⁶ Finally, as done for 3D – SALR systems using grand canonical MC simulations,¹⁰⁴ an in-depth thermodynamic characterization of the transitions between the different 2D structural patterns could be another valuable follow-up of this work.

Author contributions

GN, JKGD and CV conceived the research project. ZT developed the simulation code, conducted the simulation, and analyzed the obtained data. GN and ZT conducted the 2nd-order perturbation theory calculations. All authors participated in discussing the results and writing the manuscript.

Data availability

The simulation code and cluster analysis scripts used to generate the data in this paper are available in the referenced source.¹⁰⁵

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

ZT and GN thank Shibananda Das (GAU Göttingen), Kai Qi (2020 X-Lab Shanghai) and Roland G. Winkler (FZ Jülich) for insightful discussions, and for sharing their cluster analysis codes. The authors gratefully acknowledge computing time granted through JARA-HPC on the supercomputer JURECA at Forschungszentrum Jülich.¹⁰⁶

References

1. A. Stradner, H. Sedgwick, F. Cardinaux, W. C. Poon, S. U. Egelhaaf and P. Schurtenberger, Nature, 2004, 432, 492.
2. Y. Liu and Y. Xi, Curr. Opin. Colloid Interface Sci., 2019, 39, 123.
3. J. Ruiz-Franco and E. Zaccarelli, Annu. Rev. Condens. Matter Phys., 2021, 12, 51.
4. H. Sedgwick, S. U. Egelhaaf and W. C. K. Poon, J. Phys.: Condens. Matter, 2004, 16, S4913.
5. A. I. Campbell, V. J. Anderson, J. S. van Duijneveldt and P. Bartlett, Phys. Rev. Lett., 2005, 94, 208301.
6. C. L. Klix, C. P. Royall and H. Tanaka, Phys. Rev. Lett., 2010, 104, 165702.
7. M. Kohl, R. Capellmann, M. Laurati, S. U. Egelhaaf and M. Schmiedeberg, Nat. Commun., 2016, 7, 1.
8. N. Krishna Reddy, Z. Zhang, M. Paul Lettinga, J. K. G. Dhont and J. Vermant, J. Rheol., 2012, 56, 1153.
9. F. Cardinaux, A. Stradner, P. Schurtenberger, F. Sciortino and E. Zaccarelli, EPL, 2007, 77, 48004.
10. Y. Liu, L. Porcar, J. Chen, W.-R. Chen, P. Falus, A. Faraone, E. Fratini, K. Hong and P. Baglioni, J. Phys. Chem. B, 2011, 115, 7238.
11. F. Cardinaux, E. Zaccarelli, A. Stradner, S. Bucciarelli, B. Farago, S. U. Egelhaaf, F. Sciortino and P. Schurtenberger, J. Phys. Chem. B, 2011, 115, 7227.
12. P. D. Godfrin, S. D. Hudson, K. Hong, L. Porcar, P. Falus, N. J. Wagner and Y. Liu, Phys. Rev. Lett., 2015, 115, 228302.
13. P. D. Godfrin, P. Falus, L. Porcar, K. Hong, S. D. Hudson, N. J. Wagner and Y. Liu, Soft Matter, 2018, 14, 8570.
14. J. Riest, G. Nägele, Y. Liu, N. J. Wagner and P. D. Godfrin, J. Chem. Phys., 2018, 148, 065101.
15. E. Yearley, P. Godfrin, T. Perevozchikova, H. Zhang, P. Falus, L. Porcar, M. Nagao, J. Curtis, P. Gawande, R. Taing, I. Zarraga, N. Wagner and Y. Liu, Biophys. J., 2014, 106, 1763.
16. P. D. Godfrin, I. E. Zarraga, J. Zarzar, L. Porcar, P. Falus, N. J. Wagner and Y. Liu, J. Phys. Chem. B, 2016, 120, 278.
17. A. Chowdhury, N. Manohar, G. Guruprasad, A. T. Chen, A. Lanzaro, M. Blanco, K. P. Johnston and T. M. Truskett, J. Phys. Chem. B, 2023, 127, 1120.
18. I. Alshareedah, T. Kaur, J. Ngo, H. Seppala, L.-A. D. Kounatse, W. Wang, M. M. Moosa and P. R. Banerjee, J. Am. Chem. Soc., 2019, 141, 14593.
19. D. Pini, A. Parola and L. Reatto, J. Phys.: Condens. Matter, 2006, 18, S2305.
20. G. Foffi, G. D. McCullagh, A. Lawlor, E. Zaccarelli, K. A. Dawson, F. Sciortino, P. Tartaglia, D. Pini and G. Stell, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 031407.
21. A. L. Benavides and A. Gil-Villegas, Mol. Phys., 1999, 97, 1225.
22. G. Jiménez, S. Santillán, C. Avendaño, M. Castro and A. Gil-Villegas, Oil Gas Sci. Technol., 2008, 63, 329.
23. P. D. Godfrin, N. E. Valadez-Pérez, R. Castañeda-Priego, N. J. Wagner and Y. Liu, Soft Matter, 2014, 10, 5061.
24. A. J. Archer, D. Pini, R. Evans and L. Reatto, J. Chem. Phys., 2007, 126, 014104.
25. B. Chacko, C. Chalmers and A. J. Archer, J. Chem. Phys., 2015, 143, 244904.
26. A. J. Archer and N. B. Wilding, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 031501.
27. A. J. Archer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 031402.
28. S. Das, J. Riest, R. G. Winkler, G. Gompper, J. K. Dhont and G. Nägele, Soft Matter, 2018, 14, 92.
29. Y. Zhuang and P. Charbonneau, J. Phys. Chem. B, 2016, 120, 6178.
30. Y. Zhuang and P. Charbonneau, J. Phys. Chem. B, 2016, 120, 7775.
31. N. G. Almarza, J. Pekalski and A. Ciach, J. Chem. Phys., 2014, 140, 164708.
32. M. Litniewski and A. Ciach, J. Chem. Phys., 2019, 150, 234702.
33. J. A. Bollinger and T. M. Truskett, J. Chem. Phys., 2016, 145, 064902.
34. M. Seul and D. Andelman, Science, 1995, 267, 476.
35. D. Andelman and R. E. Rosensweig, J. Phys. Chem. B, 2009, 113, 3785.
36. H. Serna, A. D. Pozuelo, E. G. Noya and W. T. Góźdź, Soft Matter, 2021, 17, 4957.
37. G. M. Grason, J. Chem. Phys., 2016, 145, 110901.
38. J. V. Selinger, Annu. Rev. Condens. Matter Phys., 2022, 13, 49.
39. A. Ciach, J. Pekalski and W. T. Gozdz, Soft Matter, 2013, 9, 6301.
40. M. G. Noro and D. Frenkel, J. Chem. Phys., 2000, 113, 2941.
41. J. Hansen, S. U. Egelhaaf and F. Platten, Phys. Chem. Chem. Phys., 2023, 25, 3031.
42. J. Riest and G. Nägele, Soft Matter, 2015, 11, 9273.
43. J.-P. Hansen and L. Verlet, Phys. Rev., 1969, 184, 151.
44. J.-P. Hansen, Phys. Rev. A: At., Mol., Opt. Phys., 1970, 2, 221.
45. G. Nägele, A. J. Banchio, M. Kollmann and R. Pesché, Mol. Phys., 2002, 100, 2921.
46. J. Gapinski, G. Nägele and A. Patkowski, J. Chem. Phys., 2014, 141, 124505.
47. R. B. Jadrich, J. A. Bollinger, K. P. Johnston and T. M. Truskett, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 042312.
48. A. Imperio and L. Reatto, J. Phys.: Condens. Matter, 2004, 16, S3769.
49. A. Imperio and L. Reatto, J. Chem. Phys., 2006, 124, 164712.
50. J. C. Armas-Pérez, J. Quintana-H, G. A. Chapela, E. Velasco and G. Navascués, J. Chem. Phys., 2014, 140, 064503.
51. B. I. Halperin and D. R. Nelson, Phys. Rev. Lett., 1978, 41, 121.
52. D. R. Nelson and B. I. Halperin, Phys. Rev. B: Condens. Matter Mater. Phys., 1979, 19, 2457.
53. K. J. Strandburg, Rev. Mod. Phys., 1988, 60, 161.
54. K. Simons and E. Ikonen, Nature, 1997, 387, 569.
55. I. Levental, K. R. Levental and F. A. Heberle, Trends Cell Biol., 2020, 30, 341.
56. E. Sezgin, I. Levental, S. Mayor and C. Eggeling, Nat. Rev. Mol. Cell Biol., 2017, 18, 361.
57. K. S. Ramamurthi, S. Lecuyer, H. A. Stone and R. Losick, Science, 2009, 323, 1354.
58. E. Merklinger, J.-G. Schloetel, P. Weber, H. Batoulis, S. Holz, N. Karnowski, J. Finke and T. Lang, eLife, 2017, 6, e20705.
59. E. J. Dufourc, J. Chem. Biol., 2008, 1, 63.
60. O. M. Schütte, I. Mey, J. Enderlein, F. Savić, B. Geil, A. Janshoff and C. Steinem, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E6064.
61. R. E. Lamerton, A. Lightfoot, D. J. Nieves and D. M. Owen, Front. Immunol., 2021, 12, 564.
62. R. A. Dwek, T. D. Butters, F. M. Platt, T. M. Cox and T. Harder, Phil. Trans. R. Soc. B, 2003, 358, 863.
63. N. Destainville, M. Manghi and J. Cornet, Biomolecules, 2018, 8, 104.
64. V. Wasnik, N. S. Wingreen and R. Mukhopadhyay, PLoS One, 2015, 10, 1.
65. J. J. Sieber, K. I. Willig, C. Kutzner, C. Gerding-Reimers, B. Harke, G. Donnert, B. Rammner, C. Eggeling, S. W. Hell, H. Grubmüller and T. Lang, Science, 2007, 317, 1072.
66. T. Gurry, O. Kahramanohgullari and R. G. Endres, PLoS One, 2009, 4, e6148.
67. N. Destainville, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 011905.
68. N. Destainville and L. Foret, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 051403.
69. P. Charbonneau and D. R. Reichman, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 050401.
70. D. F. Schwanzer and G. Kahl, J. Phys.: Condens. Matter, 2010, 22, 415103.
71. D. F. Schwanzer, D. Coslovich and G. Kahl, J. Phys.: Condens. Matter, 2016, 28, 414015.
72. C. Bores, N. G. Almarza, E. Lomba and G. Kahl, J. Phys.: Condens. Matter, 2015, 27, 194127.
73. Y. H. Liu, L. Y. Chew and M. Y. Yu, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 066405.
74. J.-X. Chen, J.-W. Mao, S. Thakur, J.-R. Xu and F.-Y. Liu, J. Chem. Phys., 2011, 135, 094504.
75. L. Q. Costa Campos, S. W. S. Apolinario and H. Löwen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 042313.
76. Z. Tan, V. Calandrini, J. K. G. Dhont and G. Nägele, to be submitted.
77. E. Mani, W. Lechner, W. K. Kegel and P. G. Bolhuis, Soft Matter, 2014, 10, 4479.
78. F. Sciortino, S. Mossa, E. Zaccarelli and P. Tartaglia, Phys. Rev. Lett., 2004, 93, 055701.
79. J. Dhont, An Introduction to Dynamics of Colloids, Elsevier Science, 1996.
80. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 1989.
81. J. C. F. Toledano, F. Sciortino and E. Zaccarelli, Soft Matter, 2009, 5, 2390.
82. G. A. Chapela, F. del Río, A. L. Benavides and J. Alejandre, J. Chem. Phys., 2010, 133, 234107.
83. J. Torres-Arenas, L. A. Cervantes, A. L. Benavides, G. A. Chapela and F. del Río, J. Chem. Phys., 2010, 132, 034501.
84. N. E. Valadez-Pérez, A. L. Benavides, E. Schöll-Paschinger and R. Castañeda-Priego, J. Chem. Phys., 2012, 137, 084905.
85. V. M. Trejos, A. Santos and F. Gámez, J. Chem. Phys., 2018, 148, 194505.
86. A. Martinez, M. Castro, C. McCabe and A. Gil-Villegas, J. Chem. Phys., 2007, 126, 074707.
87. V. M. Trejos, M. Becerra, S. Figueroa-Gerstenmaier and A. Gil-Villegas, Mol. Phys., 2014, 112, 2330.
88. V. M. Trejos, A. Martínez and A. Gil-Villegas, Fluid Phase Equilib., 2018, 462, 153.
89. J. Bialké, T. Speck and H. Löwen, Phys. Rev. Lett., 2012, 108, 168301.
90. A. Zöttl and H. Stark, Phys. Rev. Lett., 2014, 112, 118101.
91. M. Theers, E. Westphal, K. Qi, R. G. Winkler and G. Gompper, Soft Matter, 2018, 14, 8590.
92. P. S. Ruiz, Q.-L. Lei and R. Ni, Commun. Phys., 2019, 2, 1.
93. R. J. Baxter, J. Chem. Phys., 1968, 49, 2770.
94. J. Q. Broughton, G. H. Gilmer and J. D. Weeks, Phys. Rev. B: Condens. Matter Mater. Phys., 1982, 25, 4651.
95. G. P. Hoffmann and H. Löwen, J. Phys.: Condens. Matter, 2001, 13, 9197.
96. Z. Wang, A. M. Alsayed, A. G. Yodh and Y. Han, J. Chem. Phys., 2010, 132, 154501.
97. Z. Wang, W. Qi, Y. Peng, A. M. Alsayed, Y. Chen, P. Tong and Y. Han, J. Chem. Phys., 2011, 134, 034506.
98. P. Dillmann, G. Maret and P. Keim, J. Phys.: Condens. Matter, 2012, 24, 464118.
99. J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, Oxford University Press, 4th edn, 2013.
100. F. Sciortino, P. Tartaglia and E. Zaccarelli, J. Phys. Chem. B, 2005, 109, 21942.
101. T. H. Zhang, J. Klok, R. Hans Tromp, J. Groenewold and W. K. Kegel, Soft Matter, 2012, 8, 667.
102. H. He and M. F. Thorpe, Phys. Rev. Lett., 1985, 54, 2107.
103. Z. Tan, V. Calandrini, J. K. G. Dhont, G. Nägele and R. G. Winkler, Soft Matter, 2021, 17, 7978.
104. A. P. Santos, J. Pekulski and A. Z. Panagiotopoulos, Soft Matter, 2017, 13, 8055.
105. Z. Tan, 2024 DOI:10.5281/zenodo.11694635.
106. Jülich Supercomputing Centre, J. Large-Scale Res. Facil., 2018, 4, A132.

---