Communication: Pair interaction ordering in fluids with random interactions

Multi-component systems are often encountered in studies of colloidal fluids and dispersions, which are polydisperse in particle size. They are also quite common in biology: a cell contains many thousands of different types of proteins and other macromolecules that differ in their size, shape, and in their interactions. In this work, we focus on the multiplicity of interactions between the constituents’ characteristic of biological systems (we make no attempt to model real proteins) and introduce a minimal physical model in which all the particles are different (APD) in the sense that their interaction parameters are chosen at random from a given distribution. This model differs in several important ways from past theoretical studies of multi-component systems that addressed size polydispersity of colloids and used a coarse-grained description (in terms of moments of the distribution of particle sizes) of the thermodynamic properties such as free energy and phase diagrams of these systems (see, e.g., Refs. 1–3). Conversely, in our APD model, all particles have identical size parameters and differ only in their attractive interactions with other particles. We use computer simulations to obtain insights about the temperature-dependent organization of APD fluids into non-random states in which the pair interaction parameters (PIPs) of neighboring particles become strongly correlated. Because of the huge configurational space of APD systems (there are N! non-equivalent particle permutations for each spatial configuration of N particles), we preferred to avoid a priori assumptions about equilibration and used molecular dynamics (MD) rather than equilibrium Monte Carlo methods.

The MD simulations were carried out in 2D under NVT ensemble (N = 2500 particles, density ρ^* = N/area = 0.6944) using the open source package LAMMPS.⁴ All the particles have the same mass m and size σ which are set to unity, and the simulation box is periodic in both X- and Y-axes. The particles interact via Lennard-Jones (LJ) potential

which is cut and shifted to zero at a distance r_c = 2.5σ (U_LJ(r > r_c) = 0). The pair interaction parameter ϵ_ij gives the depth of the potential well and r is the separation between a pair of particles i and j. All the physical quantities are expressed in LJ reduced units.⁵ The Langevin equations of motion of the particles are integrated with a time step of δt = 0.005τ_LJ, where τ_LJ = σ(m/ϵ)^1/2 = 1 is the LJ time (we take σ = m = 1 and define ϵ = 1 as the lower bound on ϵ_ij). The friction coefficient is fixed at ζ = 1/τ_d = 1/50τ_LJ (τ_d is the viscous damping time). We equilibrate the system at high temperature and then bring the system down to the required temperature and let it relax for about 5 × 10⁴τ_LJ to steady-state in which the ensemble-averaged (potential) energy remains constant in time. All the data reported in this work were obtained in the stationary regime.

Two different distributions P(ϵ_ij) are considered. (1) Geometric mean (GM) distribution: a particle i is assigned an interaction strength ϵ_i in the range of 1–4 taken randomly from a uniform distribution. We then adopt the Berthelot mixing ansatz⁶ and define $ϵ i j = ϵ i ϵ j$⁠. The distribution P(ϵ_ij) can be computed analytically using the methods of Ref. 7 and is shown in inset (a) in Fig. 1. It is peaked at $ϵ i j mp =2.0$ and has a mean of 〈ϵ_ij〉 = 2.42.

(2) Uniform (U) distribution: instead of assigning each particle a value ϵ_i, random ϵ_ij values are assigned to all the possible pairs of particles, i.e., N(N − 1)/2 numbers are drawn randomly from a uniform distribution in the range of 1–4 (horizontal line in inset (a) in Fig. 1). Unlike GM where the identity of a particle is determined by its value of ϵ_i, in the U system, a particle i is characterized by the set of N − 1 PIPs with all other particles (j≠i) and since the probability of generating 2 identical such sets is vanishingly small, we conclude that all particles in the system are different from each other (note that since their interaction parameters are taken from a uniform distribution, they are all identical in a statistical sense). As a reference, we use a one-component (1C) system with ϵ_ij = 2.5 (the average value in the interval 1 − 4) for all the pairs.

In Fig. 1, we present the heating curves for the average potential energy per particle as a function of T for the two systems, GM and U (we used a heating-cooling rate of dT/dt ≈ 10⁻⁵1/τ_LJ at which both systems exhibit weak hysteresis upon heating and cooling—not shown). The abrupt change of potential energy with temperature is a signature of a transition from a liquid to a solid phase, characterized by the appearance of positional and orientational (hexatic⁸) ordering. A detailed study of the low-temperature solid phase is beyond the scope of this work but preliminary results indicate that the particles form a hexagonal crystal (see Fig. S1 in the supplementary material¹⁹) that is metastable with respect to PIP ordering. This concurs with the observation that in order to prevent crystallization (which always occurs in simulations of 1C LJ systems), one has to resort to mixtures of particles of different sizes⁹ (recall that our particles have identical size). The precise nature of the liquid-solid transition in 2D is still controversial even in 1C LJ systems and may be strongly affected by finite size effects.¹⁰ In this work, we concentrate on the statistical properties of APD fluids, and therefore, we focus on the temperature range T ≥ T^*, where T^* is defined as the temperature at which the radial distribution function decays over half the system size. This yields T^* ≈ 1.0 and 1.1 for the GM and the U systems, respectively, the former being close to that of the 1C system. We verified that similar results are obtained if we identify the transition temperature with the peak of the specific heat, with the drop in mean inter-particle distance or with the increase in orientational ordering (not shown). Since the difference between the systems stems from different distributions of PIPs, we expect the energy vs. temperature curves to approach each other in the large T limit (ϵ_ij ≪ k_BT) where short-range (r < σ) repulsion dominates and to diverge from each other at lower temperatures (ϵ_ij/k_BT ∼ O(1)) where the non-universal attractive part of the potential energy plays an increasingly important role; both effects are observed in Fig. 1. Inspection of inset (b) in Fig. 1 shows that the GM system has a slightly higher energy than the 1C system, possibly because the mean of its P(ϵ_ij) distribution (2.42) is slightly lower than that of 1C system (2.5). The larger slope of the potential energy vs. temperature curve for the U system means that the specific heat $C V ∝ ∂ U ∂ T V$ and the entropy change $ΔS= ∫ T 0 T ( C V / T ) dT$ are both higher for the U system. If the energies of the U and the 1C systems were comparable, this would suggest that the liquid-solid transition should take place at lower temperature in the former system, in direct contradiction with our simulation results. This, however, is not the case, and the upward shift of the transition temperature (1.1 for U vs. 0.97 for 1C) is the consequence of the large downward shift of the energy of the U system compared to the 1C observed in Fig. 1. As will be shown in the following, this shift reflects the PIP ordering of the U fluid due to the existence of a large number of states in which pairs of particles with high values of ϵ_ij cluster together, thus reducing the mean energy of the system. Although GM fluids undergo PIP ordering as well, its extent is more limited because of the constraints imposed by the peaked form of the distribution function of the pair interaction parameters (see inset (a) in Fig. 1).

We computed the probability distribution of the total energy and found that at high temperatures, the distributions for all three systems are strictly Gaussian with width obeying the thermodynamic relation between the variance of energy fluctuations and the specific heat.¹¹ As the transition temperature is approached, deviations from Gaussian behavior are observed as the correlation length approaches the finite system size,¹² but while low energy fluctuations are enhanced and high energy ones are suppressed for the U fluids, the reverse takes place for GM and 1C fluids (see Fig. S2 in the supplementary material¹⁹). This is another manifestation of higher PIP ordering in the U system.

We proceed to take a closer look into the PIP ordering of APD fluids. In the GM fluid, there is a hierarchy in the spatial organization of particles according to ϵ_i values, i.e., particles with larger values of ϵ_i tend to form clusters at intermediate T (around T^*), while particles with small ϵ_i are preferentially localized around vacancies in the dense phase (see Fig. 2(a)). This suggests that GM fluids relax to a locally ordered state in which there are strong correlations between the interaction parameters of neighboring particles. In order to characterize the local statistical properties of both GM and U fluids (in the U system, there is no single parameter ϵ_i that characterizes the ith particle), we introduce the effective interaction parameter $ϵ i eff$ of particle i that depends not only on the intrinsic interaction parameters of the particle but also on the identity of its nearest neighbors in a particular configuration of the system,

where the sum over j goes over all the n_b surrounding particles within a cutoff radius r_c = 1.7 that corresponds to the minimum between the first and the second peaks of the radial distribution function. In Fig. 2(a), we show typical snapshots of the two APD systems taken at the same normalized distance δ = (T − T^*)/T^* ≈ 0.1 from their respective transition temperatures (large-scale particle motion characteristic of a fluid can be observed in the movies in the supplementary material¹⁹). Inspection of this figure shows that the U system is more homogeneous and has a higher mean value of $ϵ i eff$ than the GM system (for the 1C system, the distribution is a delta function at $ϵ i eff =2.5$⁠). This observation is confirmed quantitatively in Fig. 2(b) where we plotted the distributions of $ϵ i eff$ for both systems. While the U distribution is nearly Gaussian and is centered about $ϵ i eff ≃2.9$⁠, the GM distribution is much broader, with a complex shape that results from the constraints imposed by the non-uniform distribution P(ϵ_ij). Further insight into the PIP ordering of the GM fluid can be gained by considering particles with ϵ_i in a particular narrow interval of values and computing the distributions of $ϵ i eff$ for three different such intervals. We find (see inset Fig. 2(b)) that the corresponding distributions are nearly Gaussian, with width increasing with ϵ_i. Since each of the peaks is much narrower than the cumulative distribution, the GM fluid behaves as a mixture of finite number (≪N) of components, which undergoes microphase separation as T^∗ is approached (see Fig. 2(a)).

We now turn to study the temperature dependence of the distributions of $ϵ i eff$⁠. Inspection of Fig. 3(a) shows that at high temperatures, the means of the U and the GM distributions of $ϵ i eff$ approach the means of the corresponding ϵ_ij distributions, 2.5 (this happens at yet higher values of T) and 2.42, respectively. As temperature is decreased, the two distributions shift to higher values of $ϵ i eff$ and, while the width of the U distribution decreases, that of the GM distribution increases. The upward shift of $〈 ϵ i eff 〉$ is much more pronounced for the U than for the GM system, a signature of more efficient PIP ordering in the former system. Preliminary results of our Monte Carlo simulations of the crystalline phase of the U system indicate that as temperature is lowered below T^∗, the peak of the distribution continues to shift to higher values of $ϵ i eff$ and the distribution becomes asymmetric with a tail towards lower values of $ϵ i eff$ (not shown). At yet lower temperatures, the distribution is expected to narrow down and approach $ϵ i eff =4$ as T → 0 and N → ∞ (this does not happen for GM systems since particles with ϵ_i < 4 cannot have $ϵ i eff =4$⁠). Since $〈 ϵ i eff 〉$ can be thought of as new temperature-dependent energy scale, it is interesting to check whether the thermodynamics of the two multi-component systems can be reduced to that of an effective one-component system¹³ by rescaling all energies and temperatures by $〈 ϵ i eff 〉 ( T )$⁠. As shown in Fig. 3(b), this ansatz appears to work for both the U and the GM systems (in the latter case, some qualitative differences with the 1C system remain upon rescaling).

In order to gain additional insight into the nature of self-organization of APD fluids, we calculated the partial radial distribution function (PRDF) g_ij(r) that gives the probability to find a particle of type j at a distance r from particle i, for GM and U fluids at their respective transition temperatures (see Fig. 4). Here, i stands for particles with $ϵ i eff >〈 ϵ i eff 〉$ and j denotes either particles in the same (diagonal PRDF) or in the complementary range $ϵ j eff <〈 ϵ j eff 〉$ (off-diagonal PRDF). Here, $〈 ϵ i eff 〉≈2.47$ (GM) and 2.97 (U). In the U system, the only difference between the diagonal and the off-diagonal PRDFs is that the first peak is slightly higher for the former, in accord with our assertion that the PIP distribution in U fluids is spatially homogeneous. In the GM system, all the diagonal peaks are higher than the off-diagonal ones, indicating that GM fluids relax into an inhomogeneous state in which regions enriched in particles with larger values of ϵ_i possess higher long-range (hexagonal) order.

We would like to stress that the distributions in Fig. 2(b), as well as all other statistical properties we considered so far, do not depend on the particular realization of the chosen ϵ_ij distribution, i.e., they are translationally ergodic.¹⁴ In order to check whether APD fluids are ergodic in the usual sense of the word, we calculated the distribution of $ϵ i eff$ by (a) sampling the instantaneous configuration of the system (ensemble average) and averaging over many instantaneous configurations in order to reduce the noise and by (b) following the trajectories of several labeled particles over large time intervals and collecting the values of $ϵ i eff$ along these trajectories. As shown in Fig. 2(b), both methods yield identical results for the U system (in the GM system, trajectory sampling will yield different results for particles with different values of ϵ_i).

To summarize, we have studied the steady-state properties of systems of particles with random PIPs. We introduced an effective interaction parameter $ϵ i eff$ that characterizes the degree of local PIP ordering and found that this parameter increases monotonically upon cooling as the freezing transition is approached (preliminary results of Monte Carlo simulations of 2D APD crystals suggest that the increase of $ϵ i eff$ saturates upon further lowering of temperature). The constraints that define this organized state confine the system to a subspace of the available configurational space but the dimensionality of this subspace appears to be sufficiently large to allow for its efficient exploration (for discussion of diffusion in rugged high-dimensional energy and fitness landscapes, see Ref. 15), and therefore, the macroscopic behavior of APD fluids is similar to that of normal fluids. The degree of PIP ordering depends on the distribution of PIPs and is much more pronounced in U than in GM systems. This is the consequence of the fact that the choice of PIPs is much more constrained in the GM than in the U system (we choose N random numbers for GM vs. N² random numbers for U) and due to the unbiased (flat) vs. peaked forms of the PIP distributions of U and GM systems, respectively.

One may wonder whether APD systems obey standard thermodynamics since the fact that all particles are distinguishable suggests that there is a non-extensive contribution to the entropy (lnN!) associated with all possible permutations of the particles. However, in accord with general arguments about the resolution of the Gibbs paradox in classical systems of distinguishable colloidal particles,^3,16 we did not observe any violations of standard thermodynamic relations. Results on the dynamics of the systems, e.g., single particle tracking (SPT) measurements which is more relevant to make connection with SPT experiments on complex (biological) systems (see, e.g., Refs. 17 and 18) will be discussed in an upcoming publication. Other distributions of pair interaction parameters and 3D APD systems will be studied in future work.

Helpful discussions with E. Braun, D. Rapaport, A. Ermolin, D. Kessler, and E. Barkai are gratefully acknowledged. We would thank I. Parshani, R. Shalev, and M. Caspi for help with the simulations. Y.R.’s research was supported by the I-CORE Program of the Planning and Budgeting committee and the Israel Science Foundation and by the US-Israel Binational Science Foundation.