Patterns with long and short-range order in monoloyers of binary mixtures with competing interactions

Abstract

Lateral microsegregation in a monolayer of a binary mixture of particles or macromolecules is studied by molecular dynamics simulations in a generic model with the interacting potentials inspired by effective interactions in biological or soft-matter systems. In the model, the energy is minimized when like particles form small clusters, and the cross-interaction is of opposite sign. We show that the laterally microsegregated components in the dense ordered phases form alternating stripes for similar densities, or the clusters of the minority component fill the hexagonally distributed voids formed in the dense phase of the majority component. A qualitative phase diagram in the plane of densities of the two components is constructed for low temperatures. An addition of the second component significantly enlarges the temperature range of the stability of the ordered phases compared to the stability of these phases in the one-component system. At higher temperatures, the disordered phase consisting of individual particles, one-component clusters and two-component super-clusters of various sizes is stable. The product kn(k), with n(k) denoting the average number of super-clusters composed of k particles, decays exponentially with k, and the inverse decay rate depends linearly on temperature.

---

I. Introduction

Complex solvents in soft- and living matter often induce complex interactions between particles or macromolecules. A notable example is the thermodynamic Casimir potential between objects present in binary or multicomponent solvents that are close to a miscibility critical point.¹ The objects attracting the same or different components of a solvent attract or repell each other, respectively, and the range of these interactions is equal to the correlation length of the concentration fluctuations in the solvent. When the macromolecules or particles are charged, electrostatic interactions with a sign opposite to the sign of the Casimir potential appear.^1–3 When the attraction and the repulsion of any origin are both present and dominate at different separations between the particles or macromolecules, complex patterns may appear.^4–16

Interestingly, multicomponent biological membranes in homothermic living organisms are close to the miscibility critical point.^17–19 It was hypothesised that the Casimir attraction may induce aggregation of the membrane proteins, but a macroscopic separation of the proteins was not observed.^19,20 The membrane proteins are charged and according to recent discoveries, the screening length can be large in ionic solutions as concentrated as in living organisms.^21–28 When the screening length is larger than the correlation length but the strength of the screened electrostatic interactions is smaller than the strength of the Casimir potential, then the interaction between charged particles or macromolecules is attractive at short and repulsive at large distances (SALR). SALR interactions would lead to formation of aggregates of the membrane proteins of a size determined by the range of the attraction, and separated by a distance determined by the range of the repulsion. Small, well separated clusters of membrane proteins are indeed observed.

The hypothesis that like membrane inclusions interact with the SALR potential, and that the interaction between oppositely charged inclusions favouring different components of the lipid bilayer is repulsive at short and attractive at large distances was not verified yet. It inspires, however, a general question of pattern formation in monolayers of binary mixtures of particles with such interactions. We try to answer this question by considering models with different ranges and shapes of the attractive and repulsive parts of the potentials.^29–33

Although the models with the above interactions are inspired by the interactions between two types of membrane inclusions, one could design soft-matter systems that mimic such type of interactions between particles. Colloidal particles interacting with the thermodynamic Casimir potential have been already studied experimentally and theoretically.^9,34–37 Moreover, different attractive forces, such as depletion or capillary interactions can be present,^{5,14,15,38–40} and repulsion may follow from dipole–dipole interactions or soft polymeric shells.^14,16,37 Thus, one can design many soft-matter systems with competing attractive and repulsive potentials. Spontaneous patterns at different length scales may find various applications, so this problem is of practical significance as well.

Mixtures of self-assembling particles attract increasing attention,^{31–33,41–45} but a complete phase diagram has not been determined yet for any particular system. In the case of the SALR interactions between like particles and cross-interaction of the opposite sign, preliminary results have been obtained for intermediate temperatures for interactions leading to large aggregates (tens of particles) in mean-field (MF) approximation and by MC simulations.³⁰ The energetically favoured patterns on the plane of chemical potentials were determined for a potential favoring small aggregates on a triangular lattice by direct calculations and MC simulations.³² The question how the structure of the ordered phases and transitions between them depend on the shape of the interactions (determining the size of the aggregates), temperature, density and composition, however, remains to a large extent open.

By comparing the snapshots in the one-component SALR model and in the considered mixture, one can see that the structure in the dilute disordered phase in the two cases is significantly different. Dispersed small one-component clusters in the first system are replaced by much larger super-clusters made of alternating clusters of the two components in the second one.^31,33,46 However, quantitative characterization of the patterns lacking long-range order in the self-assembling mixture remains a real challenge.

In this work we focus on the low-temperature ordered phases with laterally microsegregated components in the monolayer of particles with energetically favoured alternating thin stripes of the two components. We choose interactions that favour the stripes consisting of two adjacent chains of particles of the same component (a bilayer in 2D). We perform molecular dynamics (MD) simulations, and construct a qualitative phase diagram with a topology that should be common for many mixtures with interactions favouring alternating thin stripes of the two components. In addition, we focus on the high-temperature region of the phase diagram where the disordered phase is stable, and study structure of this phase by calculating the size distribution of self-assembled clusters and super-clusters. We investigate the evolution of the structure upon increasing temperature, density and mole fraction of the particles.

In Section II, we introduce the model. In Section III, we summarize the MF results for the λ-surface separating the low- and high temperature regions in the phase space, where the phases with the long-range order are and are not expected, respectively. The simulation method is described in Section IV. Our results for the patterns formed spontaneously at low T are presented in Section V A, and for the structure of the disordered phase in Section V B. The last section contains summary and concluding remarks.

II. The model

Our model is developed for large macromolecules or colloidal particles, with diameters in the range σ ∼ 10 nm–1 μm. Interactions between such large particles differ from interactions between atoms and small molecules that are very well described by the Lennard-Jones (L-J) potential. To avoid overlap of the particles, the repulsive part of the L-J potential should be replaced by a stronger repulsion at short distances (in σ-units). Moreover, particles in statistical mechanics are modeled as hard spheres with additional interactions beyond the hard core, and to compare theory with simulations, we should assume similar interactions. For the above reasons, we assume that the particles of both components have hard or nearly hard spherical cores, and the sizes of the cores are equal. As in the previous studies,^10,29 we assume the same interactions between like particles,

u₁₁(r) = u₂₂(r) = u_hc(r) + u(r), (1)

whereas for the cross-interaction we assume

u₁₂(r) = u_hc(r) − u(r). (2)

Different shapes of the interactions were considered in the case of the one-component SALR systems, including u(r) consisting of attractive square well and repulsive ramp,¹¹ or of attractive and repulsive Yukawa potentials (double Yukawa).^13,30,37u₁₁(r) equal to the sum of the L-J and Yukawa potential^29,47 or with the repulsive part of the L-J potential replaced by r⁻ⁿ with n > 12, and the exponential repulsion at larger distances⁴⁸ were also studied. The same sequence of ordered phases was obtained in each case,⁷ therefore for the qualitative analysis of the phase behavior, only major features of the interaction potential play a significant role. This universality resembles the universal topology of the phase diagram in simple fluids, where the gas–liquid coexistence has the same qualitative features for the L-J and the hard-core square-well potentials. The optimal size of the clusters and the separation between them, however, are determined mainly by the range of the attractive and repulsive parts of the potential. From the simulations point of view, it is advantegous to study small aggregates, because many of them can be formed with a resonable number of particles, and sufficiently large number of the aggreagates is necessary for determination of their distribution in space. Importantly, clusters of membrane proteins are small, and it is not obvious if the same phase diagrams are obtained in mixtures forming large and small aggregates. For the above reasons, in this work we are interested in self-assembly into bilayers or small clusters in the absence of the second component. Self-assembly into small clusters can be induced by potentials that have different shapes, provided that a deep minimum for a relatively small distance is followed by a repulsive barrier. We verified that the potential satisfying the above requirements has in particular the form

and we choose this potential to study qualitative properties of the self-assembly in the considered binary mixture.

In the mesoscopic theory, u_hc is the hard-sphere potential, u^t_hc(r) = ∞ for r < 1 and u^t_hc(r) = 0 for r > 1 that prevents from overlapping of the particle cores. We use the superscript t for the theory. To mimic the (nearly) hard cores in the MD simulations, we assume very strong repulsion for r < 1 of the form

where the superscript s is for simulations. In the above equations, the length is in units of the particle diameter σ, and the energy is in arbitrary units . The cutoff distance in simulations is r_cut = 6.75 for all u_ij(r) with i, j = 1, 2, as in our previous works.^10,47

To compare the theoretical results obtained for u^t_hc(r) with the MD simulations obtained with u^s_hc(r) given by eqn (4), we should take into account that the minimum of u_ii(r) depends on the form of u_hc. u_ii(r) takes the minimum at r = 1, and u^t_min = u(1) = −1.815 when u_hc = u^t_hc, and u_ii(r) crosses zero for r ≈ 1.2808. For u_hc approximated by eqn (4), the minimum of u_ii(r) is u_min = u_ii(1.075) ≈ −0.6393, and u_ii(r) crosses zero for r ≈ 1.28. The large difference between the minimum of u_ii(r) in theory and simulations means that the kinetic energy and the minimum of the potential energy are equal at different temperatures in the two cases. To overcome this problem in comparing the results of theory and simulations, we assume that the appropriate energy scale is the minimum of u_ii, and introduce the dimensionless temperature T* = k_BT/|u_min|, with k_B the Boltzmann constant and u_min = u^s_min ≈ −0.6393 in the simulations and u_min = u^t_min ≈ −1.815 in the theory. The potentials u_ii for i = 1, 2 and u₁₂ are shown in Fig. 1.

Fig. 1 The interaction potential u_ij with i, j = 1, 2 between like particles, u_ii(r) (eqn (1) and (3), red line) and the cross-interaction u₁₂(r) (eqn (2) and (3), blue line). Left panel: Hard core with diameter σ = 1 (dotted line), u_hs = u^t_hs. Right panel: Nearly hard core, u_hs = u^s_hs (eqn (4)), as in the simulations. Distance is in units of σ and energy is in arbitrary units.

III. Stability of the disordered phase in mean-field approximation

In theoretical studies, the stability analysis of the disordered phase in MF gives the first approximate information about the onset of the periodic ordering. In the disordered phase, the average density of both components is position independent. This phase can be inhomogeneous when the self-assembled aggragates are mobile. It looses stability when the second functional derivative of the grand potential functional of the local densities ρ₁(r), ρ₂(r),is no longer positive definite for constant functions ρ₁(r) = [small rho, Greek, macron]₁, ρ₂(r) = [small rho, Greek, macron]₂. The internal energy for the interactions having the property u₁₁(r) = u₂₂(r) = −u₁₂(r) outside the hard cores has the formwhere c(r) = ρ₁(r) − ρ₂(r), r = |r|, k = |k| and V(r) = u(r)g(r), where g(r) is the pair distribution function. V̂(k) is the Fourier transform of V(r). In the case of hard cores, g(r) = 0 for r < 1. We neglect correlations outside the cores in MF, therefore we make the approximation g(r) = θ(r − 1). The entropic contribution in (5) is approximated by the free energy of a hard-sphere mixture in the local density approximation,where ρ(r) = ρ₁(r) + ρ₂(r) and the hard-spheres packing is described by the free-energy density f_ex in the Carnahan-Starling approximation.³⁰

The period of the self-assembled patterns should correspond to the densities ρ_i(r) that minimize the internal energy. The wave number of the energetically favored density oscillations corresponds to the minimum of V̂(k) at k = k₀. In our model k₀ ≈ 1.6558 giving the period of concentration oscillations 2π/k₀ ≈ 3.8, consistent with alternating bilayers of the particles of the two components. V̂(k) for our model is shown in Fig. 2.

Fig. 2 The interaction potential times the pair distribution function in Fourier representation, V̂(k). Red solid line: u_hs = u^t_hs (hard spheres). Blue dashed line: u_hs = u^s_hs (eqn (4), nearly hard spheres) k is in σ⁻¹ units, and V̂(k) is in arbitrary units.

The boundary of stability of the disordered phase in MF, T_λ(ρ₁, ρ₂), is called λ-surface, and was determined in ref. 30. In our temperature units defined at the end of Section II it has the form:

where T* = k_BT/|u^t_min|, c̄ = [small rho, Greek, macron]₁ − [small rho, Greek, macron]₂ and [small rho, Greek, macron] = [small rho, Greek, macron]₁ + [small rho, Greek, macron]₂ are the equilibrium values of the concentration and the density of the particles in the disordered phase, and

separates the phase space ([small rho, Greek, macron]₁, [small rho, Greek, macron]₂, T*) into the high-T region where in MF the disordered phase is stable or metastable, and the low-T region where the disordered phase is unstable. In fact local fluctuations c(r) restore stability of the disordered phase beyond MF,²⁹ and the λ-surface separates the phase space regions with dominating homogeneous and inhomogeneous distributions of the particles in instantaneous states. The ordered periodic phases can be expected on the low-T side of the λ-surface shown for our model in Fig. 3.

Fig. 3 The λ-surface of the MF boundary of stability of the disordered phase (eqn (8)). On the low-T side, the disordered phase is unstable in the MF approximation, and self-assembly into alternating aggregates of the particles of the two components takes place in typical configurations. Ordered patterns can be stable on this side of the λ-surface. Number densities ρ_i are dimensionless, .

As already noted in ref. 30, the periodic order can be present at a temperature incresing linearly with [small rho, Greek, macron] for [small rho, Greek, macron]₁ ≈ [small rho, Greek, macron]₂ (Fig. 3 and eqn (8)). For increasing difference in densities of the two components, however, the temperature at the boundary of stability of the disordered phase rapidly decreases to much lower values. This means that the ordered periodic patterns in the one component system or with a small addition of the minority component can be expected only for much lower temperatures. Note that high or low temperature refers to T* = k_BT/|u_min|, i.e. it strongly depends on the strength of the interactions.

In ref. 30, the phase diagram of a similar mixture, but with much larger period of concentration oscillations was determined under the assumption that c(r) is a superposition of sinusoidal waves with the period 2π/k₀ in different directions. Such assumption for ensemble averaged concentration can be valid for thick layers of particles and for not too low temperature.¹³ For the stripes made of two adjacent chains of the same particles, however, large deviations of the concentration profiles from the sinusoidal shape are expected in the ordered phases. For this reason, the results of ref. 30 do not apply to our model, especially at low temperature.

IV. Simulation method

We are interested in spontaneous pattern formation in monolayers of particles or macromolecules anchored in lipid bilayers or adsorbed at interfaces, where small out of plane displacements are possible. To mimic such situation in simulations, we consider particles in a rectangular box L_x × L_y × L_z with L_x = L_y = L and with periodic boundary conditions along these two directions, and L_z = 2.8. Motion of the particles along the z axis is additionally strongly restricted by the external potential that plays a role analogous to anchoring of the particles to the bilayer or interface located at the (x, y) plane,

At z = 1.07, the potential u_z(z) has a deep minimum, u_z(1.07) = −0.84 (in units). As a consequence, in the considered temperature range the equilibrium states attained by the system may be treated as two dimensional ones. Our simulations were performed applying the standard molecular dynamics method.⁴⁹ The system temperature was kept constant by scaling particle velocities. The simulation procedure is described in detail in ref. 10 and 29.

We started the initial simulations at high temperature with particles distributed randomly in the simulation box, and with N and L chosen such that the two-dimensional density ρ = N/L² was somewhat smaller than the density expected for the crystal. During the simulations, the temperature was gradually decreased until a monocrystal with the two-dimensional density ρ_s larger than N/L² was formed. The obtained monocrystal occupied a part of the simulation box, and the remaining part of the box was esentially empty.

For studies of two-phase systems we chose the final state of the simulations with N₁ = N₂ = 1800 and L = 84 (giving N/L² = 0.51), where the solid phase had the characteristic pattern with the repeated motif of two adjacent chains of particles of the first component followed by two such chains of particles of the second component, and its density was ρ = 0.76. Further simulations started from this configuration, but with the size of the simulation box increased from L = 84 to L = 180. Next the particles of the second type were added gradually to the system in the way described in ref. 10, with the temperature fixed at T* = 0.156. In this way, a few different systems were obtained with N₁ = 1800, and with N = N₁ + N₂ in the range 3840 < N < 17 600. The simulations of each system created in this way were performed at T* = 0.156 and with fixed N₁, N₂. In the additional simulations, the temperature was increased until the one-phase systems were obtained. For some systems, the simulations were carried out at the temperature up to T* = 0.469.

In order to study the structure of the disordered phase, we considered the one-component clusters in the same way as in the one-component SALR,⁴⁷ and introduced the concept of two-component super-clusters. A cluster is a compact aggregate of particles; a particle belongs to a cluster when its distance from at least one particle in the cluster is smaller than r_cl = 1.2835, where r_cl is the distance for which u_ii(r) (eqn (1)) changes sign (u_ii(r_cl) = 0). The super-cluster is a group of connected clusters of different components. Because the minimum of u₁₂(r) is very flat and the range of the interactions is large, there is no unique way of defining the connectivity between clusters of different components. We arbitrarily choose the distance r_scl = 2.15 between particles belonging to different clusters as the upper limit for formation of a bond between these clusters. We verified that different reasonable values of r_scl do not influence the results on the qualitative level. According to Fig. 1, r_scl is significantly larger than the distance corresponding to the energy minimum for u₁₂(r). On the other hand, |u₁₂(r_scl)| is only slightly lower than the absolute value of u₁₂(r) at its minimum (see Fig. 1). Cluster and super-cluster size distribution was obtained by averaging over at least 10⁵ configurations.

We compared two connectivity criteria for super-cluster formation. In the first one, one bond between particles belonging to different clusters was sufficient, and in the second one, formation of at least four bonds was neccessary to classify the two considered clusters as a part of a super-cluster. In the case of at least four-bonds formation, the bonding energy between two clusters is larger than 1.3 (in units). A disadvantage of the first criterion was formation of transient bonds during time evolution and short life-time of the super-cluster consisting of two weakly connected parts. We choose for the analysis of the structure of the disordered phase the second criterion. We verified that there was no qualitative difference between the results obtained with the two connectivity criteria. A quantitative difference, however, was significant.

V. Results

A. Low temperature phase behavior

Our aim is a construction of a qualitative phase diagram with a topology common for various mixtures self-assembling into alternating stripes or other aggregates, rather than determination of a precise phase diagram for our particular form of the interactions. Here we focus on low T, where the patterns with a long-range order can appear. To study the low temperature behavior in the simulations, we consider T* = 0.156 that is close to the λ-surface for ρ₁/ρ₂ ≪ 1 (Fig. 3), i.e. ordered phases in the one-component limit can be stable.

We simulated systems with different mole fractions and different two-dimensional (2D) densities. For presentation of the results, we choose 5 representative cases. Snapshots obtained for these representative systems are shown in Fig. 4–8 and the 2D densities in the coexisting phases are given in Table 1. In each selected system, T* = 0.156 and N₁ = 1800, whereas the total number of particles from the first to the fifth system (Fig. 4–8) takes the value N = N₁ + N₂ = 3840, 5280, 10 400, 15 200, 17 600. The average concentration and density both increase from Fig. 4–8.

Fig. 4 (a) A snapshot from MD simulations at T* = 0.156 for N₁ = 1800 and N = N₁ + N₂ = 3840 projected on the planar square with L_x = L_y = 180. The regions inside the small frames are zoomed in (b) for the gas phase and in (c) for the stripe phase.

Fig. 5 (a) A snapshot from MD simulations at T* = 0.156 for N₁ = 1800 and N₁ + N₂ = 5280 projected on the planar square with L_x = L_y = 180. The regions inside the small frames are zoomed in (b) for the gas and in (c) for the stripe phase.

Fig. 6 (a) A snapshot from MD simulations at T* = 0.156 for N₁ = 1800 and N₁ + N₂ = 10 400 projected on the planar square with L_x = L_y = 180. The regions inside the small frames are zoomed in (b) for the one-component clusters and in (c) for the stripe phase. Note the larger thickness (3 particles) of the stripes of the majority component.

Fig. 7 (a) a snapshot from MD simulations at T* = 0.156 for N₁ = 1800 and N₁ + N₂ = 15 200 projected on the planar square with L_x = L_y = 180. The regions inside the small frames are zoomed in (b) for the almost one-component clusters and in (c) for the two-component hexagonal phase.

Fig. 8 (a) A snapshot from MD simulations at T* = 0.156 for N₁ = 1800 and N₁ + N₂ = 17 600 projected on the planar square with L_x = L_y = 180. The regions inside the small frames are zoomed in (b) for the almost one-component stripe and in (c) for the two-component hexagonal phase.

Table 1 Systems with snapshots presented in Fig. 4–8. N = N₁ + N₂ and N/L² are the number of particles and the density in the quasi-2D system with area 180 × 180, where the length unit is the particle diameter. In each case N₁ = 1800. ρ_s and ρ_f are the 2D densities, and ρ¹_s and ρ¹_f are the 2D densities of the first component in the coexisting solid and fluid phases, respectively. The 2D densities in the coexisting phases were estimated by direct calculation from a central part of the crystallite, and from a large portion of the fluid phase. n_i with i = 1, 2 is the average number of bonds a particle of the i-th component forms with particles of the same component inside the crystalite of the denser phase. A pair of particles of the same component is considered as bonded when the distance between them is r ≤ r_cl = 1.2835 (see Section IV for more details)

	N	N/L^2	ρ s	ρ ^1s/ρs	ρ f	ρ ^1f/ρf	n 1	n 2
Fig. 4	3840	0.119	0.79	0.50	0.004	0.0000	4.00	4.00
Fig. 5	5280	0.163	0.76	0.47	0.050	0.0000	3.84	4.00
Fig. 6	10 400	0.321	0.84	0.4	0.210	0.0005	4.00	4.67
Fig. 7	15 200	0.444	0.75	0.26	0.380	0.0300	2.96	4.47
Fig. 8	17 600	0.543	0.80	0.26	0.470	0.0200	3.10	4.63

---

In all the shown cases, a dense two-component solid phase coexists with a less dense fluid phase where only the majority component is present (Fig. 4–6) or the minority component forms a few dispersed clusters surrounded by shells made of the majority component (Fig. 7 and 8). In Fig. 4–6, the dense phase consists of alternating stripes of the two components, and in Fig. 7 and 8, clusters of the minority component occupy hexagonally distributed vacancies in the crystal of the majority component.

We interpret all patterns with alternating parallel stripes of the two components as the same phase denoted by S. With growing difference between N₁ and N₂, first holes in the stripes of the minority component are formed (Fig. 5), and next stripes of the majority component become thicker (Fig. 6). Such evolution of the pattern leads to increased stability region of the S phase on the phase diagram (Fig. 10). The mole fraction of the first component in this phase is 0.4 ≤ N₁/N ≤ 0.6 (Table 1). In the dense phase with hexagonal or locally hexagonal symmetry, the vacancies in the majority component can be fully or partially filled with the particles of the minority component, and the mole fraction of the minority component is around 1/4.

In order to verify if the two-component hexagonal phase is stable or metastable at low T, we performed additional simulations, starting from a random distribution of the particles for ρ₁ = 0.19 and ρ₁ + ρ₂ = 0.76 corresponding to the dense phase shown in Fig. 7a and c. The resulting structure was the same as the one obtained in the procedure described in Section IV for the same ρ₁ and ρ, indicating that for such or similar values of the densities, the hexagonal rather than the S phase is stable at low T.

Let us focus on the less dense phase in Fig. 4–8. First a dilute gas of the majority component coexists with the S phase (Fig. 4). When the density increases, the particles in the gas surrounding the dense phase self-assemble into clusters (Fig. 5). These clusters tend to form a locally hexagonal pattern for still larger density (Fig. 6). We expect that the equilibrium phase at still lower temperature should contain clusters with average positions of the centers of mass forming a hexagonal lattice. Individual snapshots, however, may differ from the perfect hexagonal order, and defects may be present. When N₂ further increases with fixed N₁ = 1800, the phase coexisting with the two-component hexagonal phase consists of stripes of the majority component, as shown in Fig. 8 for N = 17 600. For this dense system, the barriers that must be overcome to reach the global minimum are very high, therefore the shown snapshot may not represent the global minimum of the free energy. Similar defects were observed in various systems self-assembling into stripes or layers.⁵⁰ In Fig. 7 and 8, micelle-like clusters with a core made of the minority-component surrounded by a shell consisting of the majority-component can be seen. They resemble structural units of the dense phase, and at lower temperatures presumably condense on the crystallite of this phase in thermal equilibrium. The interface between the two phases is smoother in thermal equilibrium at low T than in our snapshots, but we did not try to reach the final stable state, since the equilibrium shape of the monocrystal is not our goal in this study. We did not simulate still larger densities.

To quantify the observed structural changes for growing N with fixed N₁, we calculated the average number of bonds, n_i, formed by a particle of the i-th species with particles of the same species. The results are shown in Table 1 for the two-component crystalite (the denser phase). For the most symmetrical case (N = 3840, Fig. 4), a particle of each component forms on average 4 bonds with particles of the same component. When the mole fraction of the minority component decreases, n₁/n₂ decreases as well, but as long as stripes are formed, this ratio differs only slightly from 1 and takes the lowest value n₁/n₂ ≈ 0.86 for N = 10 400. n₁/n₂ takes a smaller value in this case because more bonds are formed by particles of the majority component in the central line in the thicker stripes (Fig. 6c). A significant decrease of n₁ from n₁ ≈ 4 to n₁ ≈ 3 takes place at the structural change from the alternating stripes to the hexagonal pattern, and n₁/n₂ ≈ 0.66 for N ≥ 15 200 (Fig. 7). Another quantitative characterization of the patterns in Fig. 4–8 is the probability of finding k particles in a cluster of the minority component, p₁(k). In the S phase, p₁(k) takes a maximum for the number of particles in a typical stripe that is of order of N₁ divided by the number of stripes. In contrast, in the hexagonal phase p₁(k) has the form shown in Fig. 9 for N₁ = 1800, N = 17 600 and T* = 0.156.

Fig. 9 Probability p₁(k) that a cluster of the minority component inside the two-component hexagonal phase consists of k particles, in the system with N₁ = 1800, N = 17 600 and T* = 0.156, with the snapshot shown in Fig. 8.

The simulation results and the results obtained for the triangular-lattice model,³² as well as the already established phase behavior of the one-component SALR systems^{6–8,11,13,37} suggest that the low-temperature diagram should have the qualitative shape shown in Fig. 10. The snapshots in Fig. 4–8 represent evolution of the pattern along a vertical line in the diagram in Fig. 10 for fixed ρ₁ smaller than in the one-phase regions.

Fig. 10 Schematic phase diagram at fixed low temperature. ρ₁ and ρ₂ are the densities of the first and the second component. In the white regions, single phases are stable. Dark brown regions and light patterned regions represent three-phase and two-phase equilibria, respectively. Structure of the phases is illustrated in the pictures, with red and black circles representing particles of the first and the second component, respectively.

The one-component ordered phases loose stability at relatively low T (see Fig. 3). Above their melting temperature, the dense two-component phases coexist with the disordered phase, and at this coexistence become less ordered and less dense (see Fig. 11 for T* = 0.28). The structure of the disordered phase, however, becomes more complex and self-assembly into large two-component super-clusters takes place. Further heating leads to melting of the periodic patterns made by the two components.

Fig. 11 Snapshot from MD simulations at T* = 0.28 for N₁ = 1800, N = 5280 projected on the planar square with L_x = L_y = 180.

To study the transitions to the disordered phase, we considered heat capacity and compressibility as functions of T. The heat capacity per particle, N⁻¹d〈E〉/dT, is shown in Fig. 12 with the derivative approximated by ΔE/ΔT*, with ΔT* = 0.0313 or ΔT* = 0.0156 far from or close to the maximum, respectively. For a comparison, we calculated ΔE/ΔT* for N = 3840 and two temperatures, T* = 0.156, i.e. for the two-phase system as in Fig. 4, and T* = 0.313, i.e. close to the phase transition to the disordered phase, as suggested by simulations. We obtained N⁻¹ΔE/ΔT* = 4.29 and N⁻¹ΔE/ΔT* = 34.6 for the lower and the higher temperature, respectively. The much larger value of N⁻¹ΔE/ΔT* for T* = 0.313, and the plots in Fig. 12 indicate that the heat capacity takes a maximum when the long-range order starts to disappear, causing large energy fluctuations.

Fig. 12 Heat capacity per particle for N₁ = 1800 and N = 5280, 8000, 17 600 as a function of temperature (in k_BT/|u_min| units). Lines are to guide the eye.

In Fig. 13 we show the integral of the correlation function for the number density ρ = ρ₁ + ρ₂. The density-density correlation function h(r) = g(r) − 1 can be expressed in terms of the pair distribution functions g_ij(r) between the components i, j according to g(r) = g₁₁(r)x₁² + g₂₂(r)x₂² + 2g₁₂(r)x₁x₂, where x_i is the mole fraction of the i-th component. The integral of h(r),

corresponds to the static structure factor for vanishing wavevector, S(0). χ_T in the above equation denotes the compressibility. Since we consider a quasi-two dimensional system, we used the two-dimensional integral in (11).

Fig. 13 The function G, equal to the integral of the correlation function, and related to the compressibility according to eqn (11) for N₁ = 1800 and N = 5280, 8000, 17 600 as a function of temperature (in k_BT/|u_min| units). Lines are to guide the eye.

The maxima of N⁻¹ΔE/ΔT* and G in Fig. 12 and 13 suggest a phase transition for ρ₁/ρ ≈ 0.34 and ρ₁/ρ ≈ 0.25, but for ρ₁/ρ ≈ 0.1 no pronounced maxima can be seen in the shown temperature interval. When ρ₁/ρ ≈ 0.34, the transition occurs at higher temperature than for ρ₁/ρ ≈ 0.25. These results agree with theoretical expectations of increasing temperature interval of the stability of the ordered phase with decreasing |ρ₁ − ρ₂| (see Fig. 3).

B. Structure of the disorderd phase

Characteristic snapshots in the disordered phase are shown for N₁ = 1800 and N = 8000 in Fig. 14 and 15 for T* = 0.438 and T* = 0.313, respectively. In addition to single particles of each component and small clusters of the majority component, aggregates of alternating one-component clusters of various sizes are present. An example of such an aggregate that we call a super-cluster is shown in Fig. 15. For comparison, in Fig. 16 we show a snapshot of a one-component system (N₁ = 0) with the same N and T* as in Fig. 15. One can see by visual inspection that in the one-component SALR system the size polydispersity of the aggregates is much smaller than in our mixture, and no clusters as large as the super-cluster shown in Fig. 15 are present.

Fig. 14 Snapshot from MD simulations at N₁ = 1800, N = 8000, T* = 0.438 projected on the planar square with L_x = L_y = 180.

Fig. 15 (a) a snapshot from MD simulations for N₁ = 1800, N = 8000 and T* = 0.313 projected on the planar square with L_x = L_y = 180. Note the largest super-cluster in this configuration that consists of 280 particles and is located inside the small frame. (b) The above super-cluster zoomed in. The super-clusters are identified based on the four-bond criterion (Section IV). If the one-bond criterion would be used, the particles located close to the super-cluster would be identified as its part. However, life time of such weakly connected super-clusters would be very short.

Fig. 16 Snapshot from MD simulations at T* = 0.313 for N₁ = 0, N = 8000 (one-component SALR system) projected on the planar square with L_x = L_y = 180. Total number of particles and T* are the same as in Fig. 15.

By visual inspection one can easily notice different distribution of the particles in Fig. 14 and 15. It is not easy, however, to see the difference between the cluster-size distribution of the minority and the majority component just by looking at individual snapshots. To quantify the size distribution of the clusters, we first consider each component separately and compute the number of clusters consisting of k particles, n(k), as well as kn(k), i.e. the number of particles of the considered component belonging to a cluster consisting of k particles of the same component. kn(k) is presented for N₁ = 1800, N = 8000 and two temperatures, T* = 0.438, 0.313, for each component in Fig. 17a and b. For comparison, we show in Fig. 17ckn(k) for T* = 0.313, N = 8000 and N₁ = 0. The snapshots corresponding to the histograms in Fig. 17a–c are shown in Fig. 14–16, respectively.

Fig. 17 The number of particles, kn(k), belonging to a one-component cluster consisting of k particles, divided by the number N_i of the particles of this component in the system. In each case, N = N₁ + N₂ = 8000 and (a) N₁ = 1800, T = 0.438, (b) N₁ = 1800 and T = 0.313, and (c) N₁ = 0, and T = 0.313. Red and black symbols correspond to the minority (i = 1) and the majority (i = 2) component, respectively.

At the higher T*, single particles of the minority component dominate, kn(k) decreases from the maximum at k = 1, and clusters larger than k = 6 are absent. On the other hand, particles of the majority component are most probably members of 5-particle clusters, and clusters consisting of up to 15 particles can be seen. At the lower T*, self-assembly of the minority component into 3- or 4-particle clusters is more probable than formation of dimers or monomers, kn(k) has a minimum at k = 2 and the largest (rare) clusters consist of 8 particles. While the histograms of the minority component change with decreasing T qualitatively, the histograms of the majority component remain similar, except that much fewer particles of the majority component remain in the form of monomers at the lower temperature. By comparing Fig. 17b with Fig. 17c, we can see similarity between the distribution of the particles of the majority component with the distribution of the particles in the one-component system.

Similarity between the histograms for the second component in Fig. 17b and c disagrees with the visual impression of significantly different distribution of the particles in Fig. 15 and 16. This is because the distribution of the one-component clusters does not give enough information about their connectivity, clearly seen in the snapshots.

In Fig. 18 we compare the potential energy per particle, U/N, in the mixture with N₁ = 1800 and N = 8000, with U/N in the one-component system (N₁ = 0) with the same total density. Super-cluster formation present only in the former case leads to a significant decrease of the potential energy in the whole temperature range. Because of this effect of super-clusters on the energy, the ordered phases in the mixture can be stable at temperatures higher than in the one-component system.

Fig. 18 Potential energy per particle (in |u_min| units) for N₁ = 1800 and N = 8000 (red circles), and for a one-component system (N₁ = 0) at the same total density (black squares) as a function of temperature (in k_BT/|u_min| units). Note the inflection point in the red curve at T* between 0.25 and 0.3, where the heat capacity takes a maximum, and the order–disorder transition occurs.

Unfortunately, there is no unique definition of the two-component super-cluster, because there is no unique way of defining the connectivity between the compact one-component clusters, when the minimum of the energy of a pair of different particles is very flat (Fig. 1). We choose for the analysis the presence of at least four bonds between particles belonging to different one-component clusters as the criterium that these clusters belong to a super-cluster. For the details of the super-cluster definition see Section IV. We use the same notation, n(k), for the number of super-clusters consisting of k particles.

In order to analyse the super-clusters distribution, in Fig. 19 we plot the function

P(k) = ln(kn(k)) (12)

for the same thermodynamic states as in Fig. 17a and b. Linear dependence of P(k) on k is clearly seen for suffciently large k in Fig. 19, suggesting that n(k) has the following functional formand that the probability of finding a particle in the super-cluster consisting of k particles (proportional to kn(k)) has the exponential form,

p(k) = A exp(−αk), (14)

where A is the normalization constant. A ≈ α for integer k when α ≪ 1. The exponential form means that p(k) = p(k₁)p(k − k₁)/A for any k₁ ≤ k/2, i.e. a formation of a large super-cluster occurs with a probability that is approximately a product of probabilities of findig super-clusters consisting of k₁ and k − k₁ particles, times α⁻¹ ≫ 1. Since 1/α roughly measures the size of the largest super-clusters, for k ∼ 1/α we have . Thus, the large super-clusters appear and disappear by forming and breaking bonds between smaller super-clusters almost randomly.

Fig. 19 P(k) (eqn (12)) for N₁ = 1800 and N = 8000. (a) T* = 0.438 and (b) T* = 0.313.

We should note that P(k) depends on the definition of the connectivity between the clusters quantitatively, but not qualitatively, i.e. we get the linear dependence of P(k) on k with different slopes for different connectivity criteria. Importantly, for given connectivity criterium α takes the same value for different sufficiently large values of N. We obtained the same results for α for N = 8000 and N = 32 000.

The size of the largest super-clusters is of order of 1/α that depends on the thermodynamic state. We considered a number of systems with different T*, N and N₁, and determined α from the slopes of P(k). α is shown in Fig. 20 for N₁ = 1800 and N = 5280, 8000 as a function of temperature, and in Fig. 21 for T* = 0.313 as a function of the molar fraction N₁/N for fixed number of particles N = 8010, and as a function of N/L² for the mole fraction N₁/N = 0.3. The results of simulations (symbols) perfectly fit the stright lines in Fig. 20, clearly indicationg the linear dependence of α on T*.

Fig. 20 The parameter α in eqn (14) for N = 8000 (blue symbols) and N = 5280 (red symbols) as a function of temperature T* = k_BT/|u_min|. Dashed lines are the stright lines fitted to our reults.

Fig. 21 The parameter α in eqn (14) as a function of the mole fraction N₁/N for N = 8010 (N/L² = 0.19) and T* = 0.313 (left) and as a function of N/L² for the mole fraction N₁/N = 0.3 and T* = 0.313 (right).

VI. Conclusions

We investigated patterns with long- and short-range order formed spontaneously at low and high temperature, respectively, in a monolayer of a binary mixture of particles or macromolecules. In the studied generic model, alternating stripes of the two components are energetically favoured. We did not determine the full (ρ₁, ρ₂, T*) phase diagram, but from our present and from the previously obtained results, we can draw general conclusions concerning evolution of the pattern formation.

The phase coexistences expected at very low temperature are shown schematically in Fig. 10. The ordered dense phases of microsegregated components consist of alternating stripes of the two components or of clusters of the minority component filling the hexagonally distributed voids in the crystal of the majority component. These two-component phases coexist with the phases formed by the majority component. For increasing density of the majority component, the phases coexisting with the two-component crystals have the structure predicted previously for the one-component SALR systems (Fig. 4–8). We expect such phase behavior for various mixtures with energetically favored alternating thin stripes of the two components.³² When the interactions favor much wider stripes, however, at low temperatures the microsegregation into stripes takes place, and only at intermediate temperatures the hexagonal pattern occurs.³⁰ Thus, the phase diagram depends significantly on the width of the self-assembled aggregates that in turn depens on the shape of the interactions.

The one-component ordered phases melt at relatively low temperature upon heating the system (Fig. 3), and the two-component dense phases coexist with the disordered phase. At the coexistence with the disordered phase, the phases with periodically distributed components become less dense and less ordered, and the structure of the disordered phase becomes quite complex (Fig. 11). Further heating leads to transition between the hexagonal and disordered phases, and finally between the stripe and the disordered phases (Fig. 12 and 13).

The snaphots of the disordered phase (Fig. 11, 14 and 15) show a complex pattern that looks like a rather chaotic mixture of individual particles, one-component clusters and two-component super-clusters with various sizes. The super-cluster size distribution, however, is surprisingly simple, kn(k) ∝ exp(−αk), where n(k) is the average number of super-clusters consisting of k particles. Moreover, the parameter α is proportional to temperature (Fig. 20), i.e. the size of the largest super-clusters, ∼1/α, is inversely proportional to temperature. Super-clusters of different sizes behave like different components in a multicomponent mixture, and the equilibrium between assembly and disassembly of the super-clusters resembles equilibrium in chemical reactions in such a mixture.

From the fundamental point of view, it is important that only two types of ordered patterns with laterally microsegregated components can appear in monolayers of the considered binary mixtures. While the symmetries of the 2D ordered phases in the one-component SALR and in our mixture are the same, the structure of the disordered phase in the two cases is significantly different. It is intriguing that the size distribution of the super-clusters in the considered mixture has a very simple exponential form.

From the point of view of potential applications, it is important that it should be easier to obtain spontaneously ordered patterns in this type of mixtures than in the absence of the second component, due to the much larger temperature range of stability of the phases with the long-range order.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study are available within the article and from the corresponding author upon reasonable request.

Acknowledgements

AC and ML gratefully acknowledge the financial support from the European Union Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 734276 (CONIN). WTG would like to acknowledge the support from the NCN grant no. 2018/30/Q/ST3/00434.

References

1. C. Hertlein, L. Helden, A. Gambassi, S. Dietrich and C. Bechinger, Nature, 2008, 451, 172.
2. A. Gambassi, A. Maciołek, C. Hertlein, U. Nellen, L. Helden, C. Bechinger and S. Dietrich, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 061143.
3. F. Pousaneh, A. Ciach and A. Maciołek, Soft Matter, 2012, 8, 7567.
4. A. Stradner, H. Sedgwick, F. Cardinaux, W. Poon, S. Egelhaaf and P. Schurtenberger, Nature, 2004, 432, 492.
5. A. I. Campbell, V. J. Anderson, J. S. van Duijneveldt and P. Bartlett, Phys. Rev. Lett., 2005, 94, 208301.
6. A. J. Archer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 031402.
7. A. Ciach, J. Pękalski and W. T. Góźdź, Soft Matter, 2013, 9, 6301.
8. N. G. Almarza, J. Pękalski and A. Ciach, J. Chem. Phys., 2014, 140, 164708.
9. K. Marolt and R. Roth, Phys. Rev. E, 2020, 102, 042608.
10. M. Litniewski and A. Ciach, Molecules, 2021, 26, 4532.
11. Y. Zhuang, K. Zhang and P. Charbonneau, Phys. Rev. Lett., 2016, 116, 098301.
12. B. A. Lindquist, R. B. Jadrich and T. M. Truskett, Soft Matter, 2016, 12, 2663–2667.
13. D. Pini and A. Parola, Soft Matter, 2017, 13, 9259.
14. M. Rey, M. A. Fernandez-Rodriguez, M. Karg, L. Isa and N. Vogel, Acc. Chem. Res., 2020, 53, 414.
15. V. S. Grishina, V. S. Vikhrenko and A. Ciach, J. Phys.: Condens. Matter, 2020, 22, 1215.
16. A. Rauh, M. Rey, L. Barbera, M. Zanini, M. Karg and L. Isa, Soft Matter, 2017, 13, 158.
17. S. L. Veatch, O. Soubias, S. L. Keller and K. Gawrisch, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 17650–17655.
18. B. B. Machta, S. L. Veatch and J. P. Sethna, Phys. Rev. Lett., 2012, 109, 138101.
19. T. R. Shaw, S. Ghosh and S. L. Veatch, Annu. Rev. Phys. Chem., 2021, 72, 51.
20. S. Rozovsky, Y. Kaizuka and J. T. Groves, J. Am. Chem. Soc., 2004, 127, 36.
21. A. M. Smith, A. A. Lee and S. Perkin, J. Phys. Chem. Lett., 2016, 7, 2157.
22. A. Lee, C. S. Perez-Martinez, A. M. Smith and S. Perkin, Phys. Rev. Lett., 2017, 119, 026002.
23. T. Groves and S. Perkin, Faraday Discuss., 2024, 253, 193.
24. J. Zeman, S. Kondrat and C. Holm, J. Chem. Phys., 2021, 155, 204501.
25. S. Kumar, P. Cats, M. B. Alotaibi, S. C. Ayirala, A. A. Yousef, R. van Roij, I. Siretanu and F. Mugele, J. Colloid Interface Sci., 2022, 622, 819.
26. S. Safran and P. Pincus, Soft Matter, 2023, 19, 7907.
27. A. Ciach and O. Patsahan, J. Mol. Liq., 2023, 377, 121453.
28. J. Yang, S. Kondrat, C. Lian, H. Liu, A. Schlaich and C. Holm, Phys. Rev. Lett., 2023, 131, 118201.
29. O. Patsahan, M. Litniewski and A. Ciach, Soft Matter, 2021, 17, 2883.
30. O. Patsahan, A. Meyra and A. Ciach, Soft Matter, 2024, 20, 1410–1424.
31. A. Ciach, A. de Virgiliis, A. Meyra and M. Litniewski, Molecules, 2023, 28, 1366.
32. A. De Virgiliis, A. Meyra and A. Ciach, Curr. Issues Mol. Biol., 2024, 46, 10829–10845.
33. A. de Virgiliis, A. Meyra and A. Ciach, Molecules, 2024, 29, 1512.
34. V. D. Nguyen, S. Faber, Z. Hu, G. H. Wegdam and P. Schall, Nat. Commun., 2013, 4, 1584.
35. O. A. Vasilyev, E. Marino, B. B. Kluft, P. Schall and S. Kondrat, Nanoscale, 2021, 113, 6475.
36. E. Marino, O. A. Vasilyev, B. B. Kluft, M. J. B. Stroink, S. Kondrat and P. Schall, Nanoscale Horiz., 2021, 6, 751.
37. K. Marolt, M. Zimmermann and R. Roth, Phys. Rev. E, 2019, 100, 052602.
38. R. Roth, R. Evans and S. Dietrich, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2000, 62, 5360.
39. Y. Mao, M. Cates and H. Lekkerkerker, Phys. A, 1995, 222, 10–24.
40. P. A. Kralchevsky and K. Nagayama, Adv. Colloid Interface Sci., 2000, 85, 145–192.
41. C. A. Ferreiro-Rangel and M. B. Sweatman, Mol. Phys., 2018, 116, 3231–3244.
42. J. Tan, N. D. Afify, C. A. Ferreiro-Rangel, X. Fan and M. B. Sweatman, J. Chem. Phys., 2021, 154, 074504.
43. G. Munao, S. Prestipino, J.-M. Bomont and D. Costa, J. Phys. Chem. B, 2022, 126, 2027.
44. S. Prestipino, D. Pini, D. Costa, D. Malescio and G. Munao, J. Chem. Phys., 2023, 159, 204902.
45. G. Munao, D. Costa, D. Malescio, J.-M. Bomont and S. Prestipino, Phys. Chem. Chem. Phys., 2023, 25, 16227.
46. A. Ciach, O. Patsahan and A. Meyra, Condens. Matter Phys., 2020, 23, 23601.
47. M. Litniewski and A. Ciach, J. Chem. Phys., 2019, 150, 234702.
48. A. de Candia, E. D. Gado, A. Fierro, N. Sator, M. Tarzia and A. Coniglio, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 010403(R).
49. M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids, Clarendon Press, Oxford, 1990.
50. M. Seul and D. Andelman, Science, 1995, 267, 476.

---