Cluster Morphology of Colloidal Systems With Competing Interactions

Reversible aggregation of purely short-ranged attractive colloidal particles leads to the formation of clusters with a fractal dimension that only depends on the second virial coefficient. The addition of a long-ranged repulsion to the potential modifies the way in which the particles aggregate into clusters and form intermediate range order structures, and have a strong influence on the dynamical and rheological properties of colloidal dispersions. The understanding of the effect of a long-ranged repulsive potential on the aggregation mechanisms is scientifically and technologically important for a large variety of physical, chemical and biological systems, including concentrated protein solutions. In this work, the equilibrium cluster morphology of particles interacting through a short-ranged attraction plus a long-ranged repulsion is extensively studied by means of Monte Carlo computer simulations. Our findings point out that the addition of the repulsion affects the resulting cluster morphology and allows one to have a full control on the compactness or fractal dimension of the aggregates at a given thermodynamic condition. This allows us to manipulate the reversible aggregation process and, therefore, to finely tune the resulting building blocks of materials at large length scales.


INTRODUCTION
Many body systems composed of particles interacting via a short-range attraction and a long-range repulsion (SALR) have attracted attention in the past decade due to their rich phase behavior not observed for systems interacting with purely attractive potentials. One of these is the clustered fluid phase resulting from the frustration of the gas-liquid separation. This phase has been observed in experiments, as well as in molecular simulations, see, e.g., references [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. Understanding the clustering of particles is necessary to explain changes in the dynamics and rheology of colloidal systems with competing interactions, see [15,16].
The competition between attraction and repulsion over different length scales makes possible to have different phase diagrams depending on the parameters that control both the range and strength of each potential feature. [17] presented an empirical classification for systems with competing interactions based on the range of the repulsion and attraction. When the attraction range is smaller than 20% of the particle diameter σ and a repulsion range is slightly larger than that of the attraction, these systems are considered as type I system. Particles in type II systems have a longer range repulsion while the attraction range is similar to that of the type I systems. However, the range of the repulsion is still at the order of σ. Both the type I and II SALR systems have been studied see Table 1 in reference [17]. Type I and II systems present a similar phase diagram, which includes phases as: the dispersed fluid, clustered fluid, and percolated states. One difference of these two types of systems is that type I is believed to present a gas-liquid separation below the cluster temperature if the repulsion range is short enough, as is speculated by [17]. In type II systems, the repulsion prevents the equilibrium clustering of particles beyond some limit, see, e.g. [3]. In reference [11] it is found that it is possible to find a fluidcrystal coexistence in type II systems. Finally, a type III system the attraction range is much longer than that of type I and II systems, and comparable to σ with an even longer-ranged repulsion. For those systems, some theoretical study show that there can be some very interesting gas-liquid transition behavior, see, e.g., [18,19].
In this work, we focus on the cluster morphology of type I systems, since the features of the interaction are more similar to the estimated potential for colloidal systems and protein solutions, see for example references [2,4,5,[20][21][22][23]. Besides, the structural properties of clusters have not been studied systematically yet.
The phase diagram of particles with only short-ranged attractions has been studied intensively, see e.g. [24][25][26][27][28][29][30]; and it is widely accepted that their thermodynamic behavior is determined by the extended law of corresponding states, [31]. In a similar spirit [10], proposed a generalized phase diagram for particles with competing interactions after studying type I and II systems through Monte Carlo computer simulations. According to the generalized phase diagram, particles at low and moderate concentrations form a dispersed fluid state at high temperatures, and an equilibrium clustered fluid state at temperatures below T ref c , where T ref c is the critical temperature of a reference system. Particles in the reference system interact through only the attractive part of the full SALR potential, the repulsive part, beyond σ, plays the role of a perturbation. At high enough concentrations, the dispersed fluid state percolates to form the random percolated state similar to particles with only attractions. However, when increasing the concentration, clustered fluid state can transition into a different percolating network, the cluster percolated state. Below the equilibrium clustered fluid state, it has been identified a non-equilibrium region related to the gel transition, [11]; but it is out of scope of the generalized phase diagram proposed by [10].
The manuscript is organized as follows. In Section 2, the interaction potential, as well as the computer simulation protocol are described. In Section 3, the structural signature of competing interaction systems is discussed and intimately related to the properties of clusters, which are studied in Section 5. Section 4 presents the phase diagram, where the region relevant for this work is explicitly pointed out. Finally, in section 6 the main conclusions of this work and some perspectives are presented.

COMPETING INTERACTION POTENTIAL AND MONTE CARLO COMPUTER SIMULATIONS
We have performed extensive Monte Carlo (MC) computer simulations in the canonical ensemble NVT for a system made up of spherical particles with a diameter σ, interacting through a short-ranged attraction plus a long-ranged repulsion given by the expression, the range and strength of the attraction are well determined by λ and ϵ, respectively, while the maximum repulsion strength (in units of ϵ) is given by A and the corresponding range is approximately κ −1 . The potential well is related to the reduced temperature by T * k B T/ϵ, with k B being the Boltzmann's constant and T the absolute temperature. This potential form was chosen because the contribution of the attraction, as well as the repulsion, are easily distinguishable. We have also chosen the range parameters as: λ 1.1 and κ 3 to have a simple representation of the interaction between proteins in aqueous solution, see, for example [10]. We are mainly interested in the effect of the repulsion strength on the structure and aggregate morphology. For that reason, we systematically vary the values of A 0.2, 0.5, 0.7 and 1.0. MC computer simulations were carried out at the same packing fraction, ϕ 0.1; this is related with the number density ρ * 6ϕ π . We particularly chose this particle concentration to avoid the percolation threshold, which is reached when ϕ a 0.15, see [32,33]. A typical simulation consists of 10 8 MC steps to reach the equilibrium and 10 8 extra MC steps to measure the structural properties of the system. At the lowest temperatures, additional MC steps where required to reach thermodynamic equilibrium; 10 9 MC steps before the structure was measured. In some simulations, we were not able to reach the equilibrium since particles tend to aggregate in large clusters, this situation is similar to the usual gas-liquid phase separation. In the Supplementary Material, we show the calculation of the binodal line; those results indicate that the non-equilibrium states observed in simulations are close to the binodal.
To characterize the microstructure, we calculate the so-called structure factor, S(q), which quantifies the density fluctuation correlations [34]. We have simulated a system of 11,000 particles to resolve the S(q) at small q-values, i.e., long wavelengths. Once the structure of a system was determined, the cluster size distribution, P(s), as well as the radius of gyration, R g (s), were calculated. Those quantities represent the normalized probability of finding a cluster made of s particles and the corresponding effective radius. In this second calculation, the number of particles in the simulations was decreased to 4,000, since the computational cost of larger systems is quite expensive.

STRUCTURAL SIGNATURE OF THE INTERMEDIATE RANGE ORDER STRUCTURES AND CLUSTERED PHASE
The hallmark of many SALR systems, is found in the structure factor, [10]. The S(q) presents a peak at some q-value smaller than the typical one related to the correlation between pairs of particles, April 2021 | Volume 9 | Article 637138 usually around 2π/σ. Such low-q peak has been observed in protein solutions, see, for example, [1,2,5,15,22]. However, this peak is not only due to the appearance of a clustered phase, but to a more general intermediate range order (IRO) structure, see [8][9][10]22]. For example, the IRO peak has also been observed at high concentrations where particles are already forming a percolating network, see [8]. [10] defined the clustered fluid state as that thermodynamic state of a system at which the cluster size distribution, P(s), displays a welldefined peak or maximum around s ∼ 20. This state is easily detected since P(s) decreases faster than a power law, when the system is below the percolation threshold. Also, according to [10] the clustered state is reached when the low q-peak takes a value above the critical value of 2.7. This empirical criterion works well for type II SALR systems, although there is not a physical explanation for such a value. Another criterion to identify the clustered phase has been proposed by [13] based on the value of the thermal correlation length.
In Figures 1A,B the structure factor for the systems with the lower repulsion, obtained through MC simulation (symbols) and by solving the Ornstein-Zernike equation with the Percus-Yevick closure (solid lines), is shown. In (A) it is also included the S(q) for the reference system, A 0 (dashed lines). The structure factor for A 0 and A 0.2 is equal for all temperatures considered, except for qσ < 2, that means the repulsion is too weak to affect the short-ranged correlations at separations around σ but is large enough to break up the large-ranged correlations responsible of the gas-liquid phase separation, as q goes to 0. For the system with A 0.5, the low-q peak emerges as indicated in both simulation and theory, see Figure 1B. The localization of this peak shifts to lower values as T * is lowered. It is worth to mention that the theory does not agree with the simulation at low-q values however the information provided is useful to distinguish between systems close to the phase separation, panel (A) and systems with a true clustered phase, panel (B). In panel (B) the low-q peak for the lowest temperatures is close to the limit imposed to the periodic boundary conditions, but the theory can give the trend of simulation as it does at the highest temperatures. Figures 1C,D show the cluster size distribution for systems with A 0.2 and 0.5, respectively, where P(s) is multiplied by s to avoid the biasing to monomers. The probability of finding large clusters increases as the temperature gets lower, as expected, see [8]. For example, for the system with A 0.2 at T * 0.4, that is inside the phase separation region, the distribution presents two peaks, at s 1 and ∼ 4000 (data not shown). By increasing the repulsion strength to A 0.5, P(s) indicates that the probability of finding large clusters is higher for A 0.2 than for A 0.5, since the repulsion impedes the clustering. At T * 0.375 the distribution looks different than the one obtained for other temperatures. In fact, there is a small shoulder at around 60 particles, a value larger than the reported for the clustered phase in reference [10]. Also, for this temperature, the low-q peak exceeds the value of 2.5, which is closer to the expected for a clustered fluid. Then, this case is somehow in the middle of the In (A) we also include the S(q) for the reference system the SW potential, i.e., A 0 (dashed lines). In (C,D) we also include P(s) ∼ s −2.2 (solid line), the limit case for random percolation, see [35].
Frontiers in Physics | www.frontiersin.org April 2021 | Volume 9 | Article 637138 full phase separation and the formation of a clustered phase. This observation thus indicates that for a type I SALR system, the transition temperature from a dispersed fluid state to a clustered fluid state may start shift away from the gas-liquid transition line for a reference potential system. This is qualitatively different from that of a type II SALR system. Hence, for the SALR system, the strength of the repulsion can affect the transition temperature from a dispersed fluid to the clustered fluid state. The relationship between the low-q peak and the shoulder in sP(s) is clearer for the remaining cases. Figure 2 contains the same information as in Figure 1 for systems with A 0.7 and 1.0. At low temperatures, T * < 0.375, the low-q peak of the structure factor is about 2.0, and the cluster size distribution in panels (C) and (D) present a shoulder between 20 and 30 particles, which indicates the formation of the clustered fluid. Here the maximum of the S(q) is not that close to the 2.7 threshold value, however systems are forming a clustered fluid. An important observation here is that at high s values the size distribution follows the power law P(s) ∼ s −2.2 observed in random percolation and previously noted in reference [8]. Here, clusters must form larger clusters that eventually percolate in the cluster percolated phase, see reference [10].

PHASE DIAGRAM
In Figure 3A the phase diagram for the SW fluid with an interaction range λ 1.1 is displayed; it shows the gas-liquid coexistence (blue circles) and the percolation threshold (orange diamonds). The main effect having a repulsive potential coupled with the SW is either to lower the critical temperature or to inhibit the gas-liquid phase separation, favoring the formation of the socalled intermediate range order (IRO) structure, as is discussed in references [8][9][10]. For the lowest repulsion considered here, A 0.2, we did not observed the IRO peak, as well as the clustered phase, but for the larger repulsion the simulations predict the existence of the IRO peak at all temperatures considered.
The information presented in the previous section allows us to establish the boundary between a disperse fluid, the clustered one and the phase separation. Nonetheless, the exact location of this boundary depends on the specific value of the repulsive potential, see Figure 3B. Competing interaction systems at low temperatures resemble the gas-liquid phase separation for A > 0.2, however, it is difficult to establish this from the simulation results. In the Supplementary Material, we show the calculation of the binodal corresponding for different A-values by using an approximate perturbation approach developed by [36]. This approach predicts the occurrence of the gas-liquid phase separation for all repulsive potentials considered here, which agrees well with the proposal discussed in reference [17] concerning the phase diagram of type I SALR systems. For A 0.2, the separation is predicted at T * ∼ 0.4, which agrees well with the findings of our simulations. Thus, the most interesting region to analyze the clustering is localized below the percolation threshold and at temperatures around the show the corresponding cluster size distributions. Computer simulation results are denoted by symbols, while solid lines are the theoretical solution to the Ornstein-Zernike equation with the Percus-Yevick closure. In (C,D) we also include P(s) ∼ s −2.2 (solid line), the limit case for random percolation, see [35].
Frontiers in Physics | www.frontiersin.org April 2021 | Volume 9 | Article 637138 binodal, colored region in Figure 3A. This temperature window includes states where the systems could be either a disperse or a clustered fluid and close to the expected phase separation. Figure 3C shows some snapshots of the simulated system at ϕ 0.1 and T * 0.45 for different values of A, this corresponds to the yellow star in (A). For clarity, only clusters with 10 and more particles are shown. The system with purely attractive particles, case A 0, presents the expected gas-liquid phase separation, a large cluster surrounded by small clusters and monomers. As the repulsion strength increases, the liquid phase disappears and a clustered phase emerges. It can be observed that the cluster morphology depend of the value of A as we will see below.

CLUSTER MORPHOLOGY
The cluster morphology of competing interactions systems has been studied in experiments and simulations [13]. found that equilibrium clusters are more compact than those out of equilibrium. Experiments made by [7] indicate that cluster morphology is affected by the attractive contribution. Also [37], found that equilibrium clusters tend to be elongated structures. In a previous publication, [38]; we studied the morphology of clusters made of purely attractive particles at the same ϕ-value and for temperatures above the binodal. There, we found that the morphology of clusters with more than 10 particles is determined by the strength of the attraction via the second virial coefficient, while the morphology of small clusters is insensitive to the state of the system.
In this work, we have performed a systematic study of the cluster morphology. The radius of gyration is fitted using the equation: R g Cs 1/d f , where d f is the well-known fractal dimension; this definition has been used to characterize the morphology in similar systems [39]. Figure 4 shows the results for systems at different temperatures and for repulsion strengths (A) A 0.5 and (B) 1.0. It is worth to mention that we have only considered the R g of clusters observed at least 100 times in the simulation to have a more reliable and accurate statistics. For A 0.5, systems at high temperatures follow a similar trend as the one observed for attractive particles. Small clusters have almost the same shape independent of the temperature and large clusters become more compact as the attraction becomes stronger, see Reference [38]. The changes in d f of large clusters are smaller than the ones observed in systems with purely attractive interactions and clusters keep a fractal dimension close to 2.0.
From Figure 1 it was noted the presence of small shoulder at s ∼ 60 particles, at the same value there is a small change in the FIGURE 3 | (A) phase diagram for the SW fluid with λ 1.1, our reference system. Circles correspond to the gas-liquid phase separation taken from reference [29]; diamonds point out the percolation threshold, from reference [33]; using the extended law of corresponding states to map the SW with λ 1.05 data. (B) State for particles interacting through the full SALR potential, here we identify the dispersed fluid, the clustered fluid and those states in the phase separation, see Supplementary Material. (C) Snapshots corresponding to the yellow star for different repulsion strengths. Clusters are colored randomly to be easily identified. For clarity, clusters with 10 or more particles are the only ones displayed.
Frontiers in Physics | www.frontiersin.org April 2021 | Volume 9 | Article 637138 slope of R g ; this behavior is more pronounced for the case A 1.0 where also the localization of the shoulder is shifted to s ∼ 20. Such case is shown in 4 (B), the behavior or R g is completely different from that reported in reference [38]. Firstly, the morphology of small clusters with less than approximately 15 particles becomes more compact as the temperature gets lower, the fractal dimension at the lowest temperature is around 2.5. Secondly, clusters with a larger number of particles, until ∼ 60, present a lower fractal dimension as the temperature gets lower. A similar morphology has been previously reported by [23]; they simulated a system of particles interacting through a potential fitted from experimental data for a lysozyme immersed in an aqueous solution. This means that at low temperatures compact clusters join together without merging completely, that makes possible the arrangement of these clusters in elongated structures with a lower fractal dimension, in agreement to [23]. It is important to note that for those clusters 2.5 < 2R g /σ < 4, this size is comparable with 1.5 < qσ < 2.5, which is the interval where the low-q peak was located. Finally, clusters with more than approximately 60 particles seem to have the same d f , regardless the temperature of the system. The latter could be related to the fact that at all temperatures, the S(q) is very similar at low wavelengths, i.e., q → 0, see Figure 2B. Figure 4C shows snapshots for systems with A 1.0 at different temperatures. In these cases, small clusters (5 < s ≤ 12) are colored in blue and intermediate clusters (15 ≤ s ≤ 30) in red. Small clusters at low temperatures tend to be compact structures, they are less compact as the temperature increases. The intermediate size clusters We plot, along the expected R g , the radius of gyration for clusters with a fractal dimension of 2 (blue lines), 2.5 (red lines) and 1.6 (green lines). (C) Snapshots of the system with A 1.0 at different temperatures are displayed. Clusters with 5 up to 12 particles are colored in blue and clusters with 15 ≤ s ≤ 30 in red. (D) Fractal dimension for systems with A 1.0, we classify clusters in small (s ≤ 12), intermediate (13 ≤ s ≤ 60) and large (s ≥ 100), see [40]. Vertical line is the boundary between dispersed and clustered fluid.
Frontiers in Physics | www.frontiersin.org April 2021 | Volume 9 | Article 637138 at low temperatures are still compact, but adopt an elongated configuration while at high temperatures some particles look as branches. Besides, at T * 0.4 small and intermediate clusters do not have a distinctive shape as for the lower temperatures, which agree with the fact that R g does not change its slope at 15 particles. These qualitative observations are supported by the information given in Figure 4D. There, we have computed the fractal dimension of clusters at different temperatures, they are classified as small (s ≤ 12), intermediate (13 ≤ s ≤ 60) and large (s ≥ 100), this classification agrees with the different fractal regimes shown by [40]. Since the boundary between intermediate and large clusters does not have a well-defined location, we have decided to skip some of the data. From the figure, it is clear how small clusters become more compact as the temperature gets lower than T * 0.375, i.e., the boundary between the dispersed and clustered phases. Intermediate clusters become more elongated, while the morphology of large clusters only changes slightly.

CONCLUSION
In this work, we investigated the cluster formation for systems with competing interactions. The potential parameters were chosen to represent the interactions between proteins in an aqueous solution at low salinity. The analysis of the structure factor, as well as the cluster morphology, reveals that due to the addition of the long-ranged repulsion, the gas-liquid phase separation is shifted to lower temperatures compared with that of systems with pure attractive potential. For the case A 0.2, it was evident that the phase separation happens at T * 0.4 and there is no clearer IRO peak in the structure factor. For stronger repulsion, the development of the IRO states as indicated by the IRO peak is very obvious in the structure factor. And the cluster morphology showed more complex structures due to the competition of the attraction and repulsion. At the lowest temperature analyzed here, T * 0.35, the system with A 0.5 reached a non-equilibrium state, which made difficult to state if it is a gas-liquid phase transition or a non-equilibrium cluster phase. Systems with the largest repulsion strength remained in equilibrium within the simulation window. The analysis of the cluster size percolation for systems in the clustered fluid state revealed that there is cluster formation with a preferential size in clustered fluid states. This optimal size depends also on the strength of the repulsive potential. Thus for the type I SALR systems, it is possible to produce clustered states with specific cluster size by controlling the potential parameters.
At the lowest temperatures for the stronger repulsion, particles are locally organized in very compact structures. Small clusters with less than 12 particles tend to adopt the minimum energy configuration, while intermediate size clusters have elongated configurations when these compact clusters join together with certain specific configurations. This cluster formation is different from that observed in systems made of purely attractive particles, [38]. For the SALR systems studied here, the scales at which the intermediate size clusters are observed can be related with the value at which the IRO peak appears in the structure factor. An interesting future work is the analysis of the dynamics of the clustering process at different concentrations in both the cluster fluid states and the cluster percolated states.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.