# Phase Diagram of a Reentrant Gel of Patchy Particles

###### Abstract

We study the phase diagram of a binary mixture of patchy particles which has been designed to form a reversible gel. For this we perform Monte Carlo and molecular dynamics simulations to investigate the thermodynamics of such a system and compare our numerical results with predictions based on the analytical parameter-free Wertheim theory. We explore a wide range of the temperature-density-composition space that defines the three-dimensional phase diagram of the system. As a result, we delimit the region of thermodynamic stability of the fluid. We find that for a large region of the phase diagram the Wertheim theory is able to give a quantitative description of the system. For higher densities our simulations show that the system is crystallizing into a BCC structure. Finally we study the relaxation dynamics of the system by means of the density and temperature dependences of the diffusion coefficient. We show that there exists a density range where the system passes reversibly from a gel to a fluid upon both heating and cooling, encountering neither demixing nor phase separation.

## I Introduction

In recent times we have observed a very fruitful symbiosis in the field of soft matter physics: The interplay between the synthesis of nano-and meso-sized particles with a highly versatile functionality Wang et al. (2012); Jones and Mirkin (2012); Chen et al. (2011); Ruzicka et al. (2011); Seeman (2003); Condon (2006); Stewart and McLaughlin (2004); Biffi et al. (2013); Kraft et al. (2012); Jiang and Granick (2012); Walther and Muller (2008); Chen et al. (2012); Y. Iwashita (2013) and the computational design of new model systems Glotzer and Solomon (2007); de Graaf et al. (2011); Romano and Sciortino (2011); Angioletti-Uberti et al. (2012); Smallenburg et al. (2012); Bianchi et al. (2011); Preisler et al. (2013); Hong et al. (2006); Capone et al. (2012). This combined effort can favor the emergence of new materials which have the potential for novel applications. It can also allow us to obtain a deeper understanding of some intriguing states of matter such as glasses or gels Zaccarelli (2007); Sciortino and Zaccarelli (2011). Thus the old problem of statistical physics, to connect the macroscopic properties of a system with its immutable microscopic constituents, has now been extended with the possibility of an ad hoc design of the today controllable primary particles.

This bottom-up design can capitalize on the idea of competitive interactions, i.e. the fact that the competition between distinct interaction mechanisms can cause the system to self-assemble into different local structures Sear and Gelbart (1999); Groenewold and Kegel (2001); Sciortino et al. (2004); Tlusty and Safran (2000); Russo et al. (2011). For example, the interactions between particles can be designed in such a way that at intermediate temperatures the system presents a structure favoured by entropy whereas at low temperature it has a structure that is stabilized by potential energy Angioletti-Uberti et al. (2012); Roldan-Vargas et al. (2013).

Due to their directional and selective interactions, patchy colloidal particles Jones and Mirkin (2012); Wang et al. (2012); Pawar and Kretzschmar (2008); Romano and Sciortino (2012) have been found to allow a precise control of these competitive interaction mechanisms Russo et al. (2011); Roldan-Vargas et al. (2013). Recently explored examples for this include the competition between chaining and branching in patchy colloids, where a specific design of the inter-patch interactions results in a phase diagram in which the density of the coexisting liquid approaches the density of the gas Russo et al. (2011). Another example is given by colloids that are coated with two different DNA sequences, establishing a competition between intra- and inter-particle interactions, and providing a way to tune the effective inter-particle interaction to favor crystal formation Angioletti-Uberti et al. (2012). Patchy particles can also be designed to induce, via a subtle entropy/enthalpy compensation mechanism, closed-loop phase diagrams where the low-temperature stable phase is a fluid of self-assembled weakly interacting clusters Sciortino et al. (2009, 2010); Reinhardt et al. (2011); Rovigatti et al. (2013a); Almarza (2012).

In the present work we explore one particular example of competing interactions by studying in detail the phase diagram of a binary mixture of patchy particles that has been specifically designed to create a material that continuously and reversibly passes from a fluid to a gel and from the gel to a fluid upon cooling Roldan-Vargas et al. (2013). The first species of particles, , in this binary mixture consists of particles with valence (functionality) four (i.e. particles with four attractive patches) that at intermediate temperature form an entropically favorable random tetrahedral network (the gel) structurally similar to that present in atomistic systems such as silica or silicon Binder and Kob (2011); Saika-Voivod et al. (2013). Upon a further decrease in temperature this network is fragmented by the second species () of mono-valence particles which compete for bonding with the patches of species . The possibility of fragmenting the -gel is related to the fact that if a sufficiently large number of -particles is available, each -particle can form up to four -bonds, thus leading to a energetically favorable structure. In contrast, in the fully bonded network of -particles, the number of -bonds is only two per -particle. Thus, with appropriate choices for the composition of the system and the relative strengths of two types of bonds, it can be energetically more favorable to form -bonds. As a result, the -particles progressively block the network connectivity as the temperature decreases, thus returning the system to the fluid state.

In contrast to our first study (Ref. Roldan-Vargas et al. (2013)) which focused on one fixed density and one fixed composition, here we investigate the full temperature-density-composition space of the model by Monte Carlo and event driven molecular dynamics simulations and compare these results with the ones obtained in the framework of Wertheim’s theory Wertheim (1984a, b). The use of these techniques allows us to calculate with precision the coexistence region of our system and thus to determine the range of densities and compositions at which the reversible gel is thermodynamically stable. We find that at low densities the system undergoes a phase separation into an -rich network phase and a -rich gas phase, while at high densities the system crystallizes. Finally, our study provides one of the first tests for verifying the predictions of Wertheim’s theory for binary mixtures of patchy particles, thus allowing to discuss the limits of validity of this analytical approach.

The article is organized as follows. In Section II we describe our model and provide the reader with a brief review on the relevant theoretical and methodological aspects. In Section III.A we present the study of the complete three-dimensional phase diagram of our system as obtained from Wertheim’s theory and from our grand canonical Monte Carlo simulations by using an extension of the successive umbrella sampling for binary mixtures Rovigatti et al. (2013b). Section III.B concentrates on the study of the two-dimensional temperature-density cuts of the phase diagram corresponding to the stoichiometric molar fraction, i.e. the molar composition for which there are four -particles for each -particle. In Section III.C we document the spontaneous crystallization of our system — at the stoichiometric molar fraction — in molecular dynamics simulations at high densities and intermediate temperatures. Also for this particular composition, in Section III.D we focus our attention on the dynamics through the study of the temperature and density dependences of the diffusion coefficient of the -particles. Finally, in Section IV we conclude with a summary of our main findings.

## Ii Model and Theoretical Background

### ii.1 Interaction Potential: Kern-Frenkel Model

The particles investigated in this paper are given by the well-known Kern-Frenkel model Kern and Frenkel (2003), i.e. the potential between particles results from the sum of a hard-sphere potential and an attractive directional interaction . The hard-core repulsion between two particles and is given by:

(1) |

where , with Boltzmann’s constant and the temperature, is the center-to-center distance between the particles, and denotes the contact distance between the particles, with the hard-sphere (HS) diameter of particle . The site-specific attraction between the particles is determined by the circular patches on the surface of each particle, which interact such that two particles form a bond with interaction energy if i) their centers of mass are within a maximum interaction range , and ii) the center-to-center vector between the particles passes through a patch on the surface of both particles Kern and Frenkel (2003). The size of the patches is determined by an opening angle (see Fig. 1). The attractive part of the potential energy of two particles and is thus given by

(2) |

where is a square-well attraction, given by:

(3) |

The function is defined as

(4) |

where , is a set of normalized vectors pointing from the center of particle towards the center of each of its patches, and is a unit vector in the direction of .

### ii.2 Numerical Details of the Simulation

In the binary mixture of patchy particles studied here, the first species, , has a HS diameter with four identical patches on its surface that are arranged in a tetrahedral geometry. The second species, , has a HS diameter and only one patch on its surface. The attractive patch-patch interactions are modelled by the Kern-Frenkel potential discussed above. The patch interaction energies, Eq. (3), for the - and -bonds are and , respectively. Bonds between two -particles do not occur (i.e. ). The interaction ranges for the patch-patch interaction are and and the angular width of the patch is given by and (see Fig. 1 and Ref. Roldan-Vargas et al. (2013)). Due to the chosen geometry of the patches, each one can be involved only in a single bond. In addition, the ratio has been chosen such that the contribution of the -particles to the total packing fraction, close to the stoichiometric composition, was significantly smaller than that of the -particles, while maintaining the ability of the -particles to block the formation of -bonds Barcenas et al. (2007). The values of the patchy interaction energies,, were chosen to generate at low a ground state in which each -particle is bonded to four -particles, completely saturating its patches. From now on, this low preferred state will be called flower state. The values chosen for the parameters that define the geometry of the bonding volume (i.e. and , see Appendix) reflect the requirement of having an entropy associated to the -bonding that is larger than the one corresponding to an -bond, thus favouring at intermediate temperature the formation of a highly connected network of tetrahedrally coordinated -particles.

We perform event-driven molecular dynamics (EDMD) simulations Rapaport (2009); de la Peña et al. (2007); Smallenburg and Sciortino (2013) for systems of and particles (i.e. molar composition ), for different total number densities , and , covering a wide range of the temperature for each number density. The implementation of the EDMD simulation relies on the numerical prediction of bond formation and bond breaking events, using the numerical techniques described in Ref. Smallenburg and Sciortino (2013). In this work, the mass of each particle is taken to be the same, setting the time unit of the simulation . Similarly, the moments of inertia tensors of all particles are chosen to be the same: . During the equilibration of the simulations, the temperature is controlled by an Andersen thermostat: periodically, randomly selected particles are given a new translational and angular velocity, drawn from a Maxwell-Boltzmann distribution. While measuring the dynamic properties of the system (here the diffusion coefficient) no thermostat is used, so that the total energy of the system is constant. To reduce equilibration times at low , we used standard Monte Carlo simulations () to generate initial equilibrium configurations.

To calculate numerically the phase diagram of the binary mixture we use the extension of the successive umbrella sampling (SUS) method Virnau and Müller (2004) recently developed Rovigatti et al. (2013b). The SUS method provides an efficient way to sample in a grand-canonical ensemble (at fixed chemical potential, temperature and volume) the probability that particles are present in the simulated volume by dividing the interval (where is the largest sampled density) into intervals (in which the number of particles can fluctuate only by one unit). In the extension to binary mixtures, the intervals and are partitioned into overlapping two-by-two windows ( ) and in each of these windows a grand canonical Monte Carlo (GCMC) simulation is performed at fixed and , rejecting all moves which would bring the number of particles outside the selected range. We have selected and , corresponding to a total of 24000 distinct simulations. To improve the method even further, the chemical potential of the two species and is optimized in each window. The and values are selected in a preliminary run in such a way that the states with and particles (as well as the states with and particles) are sampled with comparable probabilities. With this choice, data in all windows have approximatively the same statistics. At the end of the simulations, histogram re-weighting is performed on the probabilities of the different windows to generate a set of data corresponding to the same values of and . A splicing procedure, described in detail in Ref. Rovigatti et al. (2013b), is then used to evaluate the probability distribution function . The resulting is then analysed (after re-weighting to selected and values) to determine the density and the relative fraction of the two coexisting states Rovigatti et al. (2013b). The significant computational investment required is compensated by the possibility of calculating the full behavior in the plane at the chosen . We have repeated the procedure for three different ’s (0.09, 0.12, 0.135). Smaller ’s require longer equilibration times, prohibiting us from exploring this region.

### ii.3 Wertheim Theory

In this work, the analytical approach to the thermodynamics of our system is based on Wertheim’s thermodynamic perturbation theory which allows us to obtain an analytical expression for the Helmholtz free energy, , of pure patchy particles fluids and fluid mixtures (a detailed description can be found in Refs. Wertheim (1984a, b)). Accordingly, the total free energy per particle, , can be split into two contributions:

(5) |

where is the contribution to the total free energy due to the hard-sphere interaction whereas contains the contribution due to the attractive (bonding) patchy interaction. can be written as the sum of an ideal gas contribution, , and an excess term, , which accounts for the excluded volume due to the finite size of the particles. In case of a binary mixture with molar compositions (), the ideal term is given by

(6) |

where is the thermal volume, whose constant contribution (that here we take as ) is irrelevant for the phase behavior. Concerning the excess term, , we have used the analytical formalism developed by Mansoori and co-workers for binary mixtures of HS Mansoori et al. (1971), which takes into account the different HS diameters of the two species (the explicit form of is reported in the Appendix).

The bonding free energy, , specialized to our system, is given by

(7) |

where is the probability that a patch of type is not bonded. These probabilities, or more conveniently the probabilities that an -site is bonded can be calculated through the law of mass action de las Heras et al. (2012) (see Appendix).

To obtain the equilibrium properties of our binary system we follow Ref. de las Heras et al. (2012) and minimized the Gibbs free energy per particle, , for a fixed pressure , , and composition with respect to the total number density, . Coexisting points are then determined by imposing chemical equilibrium through the equality of the chemical potentials of both species in the two coexisting phases (thermal and mechanical equilibrium are already assured by imposing a constant and ). This condition for chemical equilibrium is geometrically implemented via a common-tangent construction of the function for the two coexisting compositions at fixed and (see Appendix).

## Iii Results

### iii.1 Phase Diagram: Composition Versus Density

In this section we first discuss the composition-density profiles of the coexistence regions at different ’s as predicted by Wertheim’s theory. In general, the coexistence region of a binary mixture defines a three-dimensional volume in the space that we present here by two-dimensional cuts in the plane for different values of .

Figure 2 shows the prediction of Wertheim’s theory for our system at different ’s. According to the theory, for the system phase separates into two phases, differing in their density and/or composition, whereas for no phase separation is predicted. The figure also shows the tie lines (see Appendix), i.e. the set of points that separate into the same two coexisting states. Vertical tie lines indicate those points of the phase diagram with a prevalent demixing, that is, the two coexisting phases only differ in their composition but not in their total number density. Horizontal tie lines indicate gas-liquid type separation, that is, the two coexisting phases differ in density but not in composition. As Fig. 2 shows, in our system the gas-liquid type separation only occurs close to (the monodisperse case).

As we see in Fig. 2 we expect at all ’s phase separation at low densities (typically ), that would essentially consist of two phases (see coexisting points), one rich in -particles (high ) whereas the other would mainly consist of -particles (). This phase separation is analogous to the one that occurs in the one component system of particles with four-patches Rovigatti and Sciortino (2011), which cluster together in order to form an -network.

At high (Figs. 2a) and b)), the substantial difference between the bonding volume associated to the -bonds, , (see Appendix) and that corresponding to the -bonds, , () plays an essential role, making the formation of -bonds entropically unfavorable as compared to the -bonds. Thus, the addition of -particles to the system only increases the total number density of the two coexisting phases, without preventing phase separation. Due to excluded-volume effects, these -particles will preferentially go to the phase where the density of -particles is low, resulting in a demixing of the two species.

The formation of clusters of -particles at these temperatures () is clearly reflected by a high probability of forming -bonds, . When this probability is given by (see Appendix and Ref. Roldan-Vargas et al. (2013)). Figures 3a) and b) show for at different densities as obtained by Wertheim’s theory and from our EDMD simulations respectively. Also shown in these graphs is the probability that a patch of type of a particle of the species is specifically bonded to a patch of a -particle. Indeed, around , shows a maximum (Fig. 3a)), indicating a high degree of percolation of the -particles. The position of this maximum, , slightly increases upon increasing density, as can be recognized from these figures.

We see that with decreasing , Figs. 2c)-f), the range in density in which phase separation is observed shrinks rapidly if is small. This behavior is related to the fact that for the system becomes dominated by a state in which for energetic reasons the -bonds occur very frequently. The emergence of this new state can be recognized from the quick increase of and the decrease of , see Fig. 3. Thus for small values, typically , the high concentration of -particles limits the clustering of -particles observed at intermediate temperatures, thus driving the system back into the homogeneous phase, Figs. 2d)-f). However, for high values of , but still low densities, the number of -particles is too small to block all of the -patches. As a result, the system will still phase separate into an -rich network phase coexisting with a phase rich in -particles that are completely blocked by four -particles, i.e. that form the local structure that we referred to as flowers (see section II.B). In this case we have a phase separation between two phases that differ significantly in composition and density. In fact, this is the only phase separation that survives at very low for and , where a high fluid (essentially a pure fully bonded network) will coexist with a low density fluid with (essentially a gas of flowers) (Figs. 2e) and f)).

Next, we compare the prediction of the Wertheim theory with the corresponding composition-density cuts through the phase diagrams as obtained by the SUS method (see section II.B). Figures 4a) and b) show this comparison for two different ’s ( and ). In both cases, the boundary of the coexistence regions as obtained from the theory is in very good agreement with the one obtained from the simulation, although at high densities some differences can be seen. (Typically for and more noticable at .) For , the lowest for which we were able to equilibrate the SUS simulations, the agreement between the boundary of the coexistence region predicted by the theory and the one by the simulation is certainly excellent, see Fig. 4a). However, even at this we see that this agreement does not hold for the tie lines connecting coexisting points at intermediate densities (). As mentioned above, these deviations occur for those states for which we have a well formed -network with a high degree of percolation, i.e. geometries for which the theory is less accurate due the presence of closed loops of particles not considered in the theory de las Heras et al. (2011). For , where the degree of percolation is slightly less pronounced, the agreement between Wertheim and SUS tie lines is recovered for the whole set of densities and compositions computed by the SUS method.

### iii.2 The Phase Diagram: Temperature Versus Density

In this section we discuss the 2D cut (temperature versus density) through the phase diagram corresponding to a fixed molar composition . This is the so-called stoichiometric molar fraction since for this value of all the -patches can be bonded to -patches, i.e. the system ground state is a pure fluid of flowers. However, we mention that this is certainly not the only interesting molar composition, since, e.g., binary mixtures of patchy particles with non-stoichiometric ratios Smallenburg et al. (2013) are a valuable model for describing vitrimers, a recently invented malleable network plastic with controlled healing properties Montarnal et al. (2011).

Figure 5 shows the coexistence region for in the -plane as predicted by Wertheim’s theory. Also included are the numerical results for the three temperatures using the SUS method. According to Wertheim’s theory, no phase separation can be expected for densities , independently of the temperature. Moreover, for all densities the system is predicted to be homogeneous if . As mentioned in the previous section, the phase diagram shows a re-entrant behavior typical of systems with competing interactions Russo et al. (2011); de las Heras et al. (2011); Angioletti-Uberti et al. (2012). At intermediate ’s () the system phase separates into a network phase mostly composed of -bonds and another phase rich in -particles. As discussed above, this separation is essentially a demixing where the two phases have a similar total number density but significantly different composition, and therefore local packing fraction. (As an example, see the states inside the coexistence region with in Figs. 2 and 4). Indeed, one should realize that the two branches limiting the coexistence region (black dots joined by solid lines in Fig. 5) in this cut at constant do not correspond to two sides of the same coexisting points since in general the state points coexisting with those that appear in Fig. 5 will have a different composition, (see Fig. 2). Therefore, the phase separation present in our system will differ from that occurring in pure systems of particles with valence four for which we will have a gas-liquid type separation. However, due to the neutral role of the -particles at high and intermediate before entering into the coexistence region, the critical points as predicted by Whertheim’s theory for our binary mixture and for a pure -system () are similar. In our case, and according to Wertheim’s theory, the high border of the coexistence region is placed around with a critical partial density , whereas for the pure system we obtained with a critical density . We mention that these latter values differ from those obtained by numerical simulation for pure systems with similar interaction parameters Foffi and Sciortino (2007), where the critical point is found at but with a significantly higher density, . It is indeed well know that Wertheim theory underestimates the density of the coexisting liquid branch Bianchi et al. (2006a).

At sufficiently low , the -interactions starts to become dominant, hence favouring the formation of inert flower structures, and thus return the system to the homogeneous state for . As already shown in the previous section, there are discrepancies between Wertheim’s theory and the SUS method at intermediate ’s () where the theory underestimates the coexistence density of the network phase. Indeed, at and the numerical simulation shows that the range of the coexistence region extents up to . On the other hand theory and simulation agree again for , i.e. when the system becomes less percolated due to the emergence of the fluid of flowers at low ().

Also shown in Fig. 5 are the percolation lines as obtained by the Flory-Stockmayer theory Flory (1953), with bonding probabilities evaluated via Wertheim’s theory (always for ). Since in our system percolation is associated to the formation of the -network (-particles merely act as blockers), the percolation temperature for a given density will be that at which , that is, our criterion to establish the percolation temperature of the AA-network is formally equivalent to that one that applies for one-component systems with four valence particles Flory (1941, 1953). We have also computed the percolation points by directly analyzing the cluster distribution from our EDMD simulations. In this case, we consider a system to be percolated if the largest cluster in the system spans the simulation box. This procedure is performed for an ensemble given by several independent configurations. The percolation locus is defined as the set of state points at which at least fifty per cent of the configurations percolate. Interestingly, the locus of the high percolation line obtained by Wertheim’s theory overestimates the temperature at which the percolation threshold is reached with respect to that obtained by our EDMD simulations (see Fig. 5, blue diamonds and orange squares, respectively). This difference is likely the result of loop-formation in the fluid (present in simulations but not in the theory). Another possible source for the deviation is the analytical overestimation of the bonding volume that we performed by using the contact value of the partial radial distribution function between -particles for the whole interaction range, (see Appendix). Interestingly, the low percolation lines as obtained from Wertheim’s theory and EDMD show a very good agreement.

a) | b) |

c) | d) |

Finally, we show in Fig. 6 four snapshots corresponding to different equilibrium configurations obtained from our EDMD simulations for and , the temperature at which the coexistence region extends to its maximum density according to Wertheim’s theory (see Fig. 5). In Fig. 6 we see a clear phase separation for densities and (Figs. 6a and b, respectively). As discussed in the previous section, the separation results in two phases, one rich in unbonded -particles whereas the other mainly consists of an -particle network. For , that is, the density at which Wertheim’s theory predicts the edge of the boundary of the coexistence region, the phase separation is still obvious (see Fig. 6c), confirming the discrepancy between the prediction of the simulation and Wertheim’s theory. Finally, for , a density in the vicinity of the edge of the boundary of the coexistence region as predicted by the SUS method, we only see a small amount of phase separation, consistent with the SUS prediction. For , the system is homogeneous for all temperatures.

### iii.3 High Density Crystal

a) | b) |
---|---|

As we have already mentioned, for , intermediate temperatures , and rather high densities our system consists of a -particle bonded network with the free smaller -particles homogeneously diffusing through the network voids (see Fig. 6d). An example of this situation is the reversible gel studied in our previous work at a constant density , where a maximally connected amorphous -network at can be reversibly melted upon both cooling or heating the system Roldan-Vargas et al. (2013). If the density is even higher (e.g. , the highest density investigated) we observe a spontaneous crystallization, at temperatures , , and .

The occurring crystallization process can be followed in simulations up to the point where most of the -particles are incorporated into a crystalline environment. The system forms a coexistence of two crystal phases: a body-centered cubic (BCC) and diamond cubic (DC) crystal phase, both known to be stable for monodisperse particles with a tetrahedral patch geometry Romano et al. (2010). The BCC phase simply consists of two interpenetrating diamond lattices in which one of these lattices stretches through the entire simulation box, while the BCC region is formed in the region where the second diamond lattice is present as well. Figure 7 shows snapshots of the system, with the two separate lattices indicated by different colors. Interestingly, only the cubic form of the diamond crystal structure was formed during the nucleation of the crystal. In all cases investigated, a cluster of the BCC crystal structure was seen to form first. Subsequently one of the two diamond lattices inside the crystal structure grows and eventually fills most of the simulation box. For temperatures between , initializing a simulation with a BCC cluster present in the box (i.e. from a configuration taken from the nucleation process at ), still results in a fully crystallized system, suggesting that (barring finite-size effects) a crystalline state is stable there. At lower ’s (), the crystal melts, and we only observe the amorphous gel structure.

The absence of the hexagonal diamond phase, which is known to have a very similar free energy Romano et al. (2011), might be explained by the presence of the BCC crystal region, which is not compatible with the hexagonal diamond crystal structure, or by the depletion effect created by the small -particles Kumar and Panagiotopoulos (2013). Indeed, the BCC part of the crystal structure always nucleates first in our simulations. The BCC crystal region typically spans a significant part of the simulation box along at least one axis, and therefore can affect the crystal structure in the entire volume. In much larger systems, stacking faults resulting in a mix between the cubic and hexagonal diamond structures might still occur far away from the BCC cluster.

It is worth noting that recent studies have shown that tetrahedral patchy particles spontaneously crystallize in a DC or DH structure when the patch opening, , is smaller than approximately . For the Kern-Frenkel model, this corresponds to . Larger opening angles stabilize the formation of glasses Romano et al. (2011), while at even larger angles the liquid becomes thermodynamically more stable than the crystal Smallenburg and Sciortino (2013). It is thus intriguing to observe spontaneous crystallization in this model (for which ) in the presence of a second component (-particles).

The emergence of the crystalline regime and its subsequent melting is also nicely captured through the temperature evolution of the partial radial distribution function of the -particles, . Figure 8 shows this evolution, where has been computed at different temperatures from our EDMD simulations data. At we vaguely see the emergence of the tetrahedral amorphous structure through an incipient second peak located around . At we clearly see the signature of the crystalline phase through the well-developed maxima and minima that extent up to large and that correspond to the different distances that define the crystal structure. At the amorphous tetrahedral gel is again recovered, with only a modest second tetrahedral peak surviving (). Finally, at (the lowest temperature investigated), we see a new peak located at around () indicating the emergence of a fluid composed of flowers Roldan-Vargas et al. (2013).

### iii.4 : Temperature and Density Dependences of the Diffusion Coefficient

Finally we discuss the relaxation dynamics of our system and its density and temperature dependence for the molar composition as obtained from our EDMD simulations. For this, we calculate the diffusion coefficient of the -particles from the long-time diffusive regime of their mean square displacement (MSD) via the Einstein relation for those densities and temperatures at which the system is not phase separated. We recall that the center of mass of a single species (here the -species) may have a non-zero velocity, which is compensated by the center of mass (CM) motion of the other species. When particles form a network, like in the case, this finite-size effect contributes to the MSD of the single species. To correct this effect we subtract the CM drift of the species in question before evaluating the MSD. In addition, to discard the trivial trend originated from the -dependence of the thermal velocity we divide by a reference diffusion coefficient (see section II.B).

Figure 9 shows the -evolution of for all the densities investigated. For , there exists a region of intermediate temperatures where the system phase separates, preventing the possibility of estimating the diffusion coefficient. The -range of the phase separation as seen by the omitted values becomes narrower upon increasing , as was already seen in the phase diagram of Fig. 5. Similarly, at large densities (), the system crystallizes at intermediate , again preventing the possibility of covering the complete -range of the diffusion coefficient. Fortunately, at intermediate densities (), it is possible to follow the entire dynamical process and its non-monotonic evolution. Figure 9 clearly shows that, on cooling, the diffusion coefficient of the -particles first decreases by four orders of magnitude and then increases going back to typical fluid-like values. Reference Roldan-Vargas et al. (2013) discusses in detail this process, associated to the progressive formation of the -network which is then replaced by the gel-fractioning induced by the progressive formation of -bonds. This behavior can be rationalized by a competitive interaction mechanism that in the present case leads to an entropically favorable -bonding at intermediate temperatures and to an energetically favorable -bonding at low temperatures.

Finally, Fig. 10 shows iso-diffusivity points Foffi et al. (2002); Zaccarelli et al. (2002) corresponding to different values. At intermediate (typically ) we find the expected behavior of a fluid-like system: The curves have a positive slope, i.e. states with the same value of the diffusion coefficient have a higher temperature if density is increased (points connected by solid lines in Fig. 10). However, at low the opposite behavior appears, that is, states sharing a common value of the diffusion coefficient have a slightly lower value as density increases (points connected by dashed lines in Fig. 10). This a priori odd behavior can be understood by considering the density dependence of (or ) at low (Figs. 3a and b). From these figures we recognize that at a fixed (low) , increases as density increases, indicating that for a fixed value of the temperature the system will be more fragmented at higher densities. This increasing degree of fragmentation causes the system to diffuse faster as increases.

## Iv Conclusions

In this work we have presented an extensive study of the phase diagram of a binary mixture of and patchy particles which was specially conceived to form a reversible gel upon both cooling and heating Roldan-Vargas et al. (2013). The design of the interaction geometry and the choice of the interaction energies results in a competing mechanism between an entropically favorable network at intermediate temperatures which is formed by -particles (the gel) and an energetically favorable fluid composed of -particles completely bonded to -particles at low temperatures. To explore the phase diagram of the system (and therefore its coexistence region), we have covered both analytically and by means of computer simulations a wide range of temperature, density, and composition.

In the analytical approach to the thermodynamics of the system we have used a hard-sphere expression for the free energy (in which we have taken into account the different hard-sphere diameters of the two species) and added a bonding free energy contribution based on Wertheim’s first-order perturbation theory Mansoori et al. (1971); Wertheim (1984a, b). This approach provided us with an analytical expression for the Helmholtz free energy and therefore with the full thermodynamics of our system. In addition, we have performed two kinds of computer simulations for the same model. On one hand, we have used grand canonical Monte Carlo simulations to numerically calculate the phase diagram of our model at different temperatures by using an extension of the successive umbrella sampling (SUS) technique for binary mixtures Rovigatti et al. (2013b). On the other hand, to study also the dynamics of our system, we have performed event-driven molecular dynamics simulations for a fixed molar composition , which represents the stoichiometric molar fraction in our system for which all the -patches can be bonded to -particles.

Our calculations have allowed us to analyse the 2D cuts, composition versus density, of the phase diagram at different fixed temperatures as obtained from Wertheim’s theory and SUS. States at which the system is phase separated usually consist of two phases that are demixed: A network phase rich in -particles and a fluid rich in -particles. For small values (), this phase separation extends to higher densities upon increasing temperature. At low , phase separation only takes place at and low densities. The coexistence regions as predicted by Wertheim’s theory and SUS showed in general a good quantitative agreement, with some deviations for those states at which we have a high degree of connectivity of the -particles. In particular, while Wertheim’s theory predicts no phase coexistence for at any , the SUS method extends the coexistence region to .

Using EDMD simulations, we have specifically studied the composition , confirming the extent of the coexistence region predicted by SUS until . In addition, we have also established the temperature above which, independent of density, no phase separation is expected, finding it to be around . Since in our binary system phase separation is mainly due to demixing, in general the boundary of its coexistence region differs from that corresponding to a pure system for which we can only have a gas-liquid type separation Foffi and Sciortino (2007); Bianchi et al. (2006a); Bianchi et al. (2008). Also for this composition, we have shown that the percolation lines as obtained from Wertheim’s theory and from our EDMD simulations are in good qualitative agreement, showing a large density range where the systems only percolates in the intermediate temperature regime. In other words, upon cooling, the system first forms a percolating network, that then fragments again upon further cooling.

For we have also investigated the spontaneous crystallization of our system at high densities () and intermediate temperatures () as well as its subsequent melting recovering the gel state upon decreasing . In this crystalline regime, the system spontaneously forms a coexistence of two crystal phases: a body-centered cubic (BCC) and a diamond cubic (DC) crystal phase. Interestingly, it appears that stacking faults associated to the hexagonal and cubic forms of diamond are not observed Romano et al. (2011).

Finally, we have studied the relaxation dynamics of our system through the diffusion coefficient of the -particles for at those temperatures and densities for which no phase separation or crystallization is present. We have shown that for densities around we can follow the complete non-monotonic temperature evolution of the diffusion coefficient.

In summary, we have presented an extensive study of the thermodynamics and dynamics of a reversible gel of patchy particles by using different numerical simulations complemented with an analytical approach. With these techniques we have explored the space that defines the three-dimensional phase diagram of our system, locating its coexistence region and describing its different phases. In particular, we have presented a promising result for future experimental realizations Stewart and McLaughlin (2004); Biffi et al. (2013); Wang et al. (2012) by demonstrating that there exists a broad density range where our system exhibits the phenomenology of a thermo-reversible gel that can be fluidized by both cooling and heating in the absence of phase separation.

## V Acknowledgements

We acknowledge support from COMPLOIDS, ERC-PATCHYCOLLOIDS-226207 and MIUR-PRIN. W. Kob acknowledges support from the Institut Universitaire de France. We thank L. Rovigatti for providing us with the code to perform and analyze the SUS simulations.

## Vi Appendix: Free Energy of a Binary Mixture of Hard-Spheres With Patchy Interactions.

As mentioned in section II.C, the Helmholtz free energy of a binary mixture of HS, , can be separated into an ideal gas contribution, , and an excess term, (see Eqs. (5) and (6)). In our study we have considered the different particle HS diameters following the approach of Mansoori and co-workers Mansoori et al. (1971) in which is written as:

(8) |

where:

(9) |

(10) |

(11) |

(12) |

and

(13) |

Here we have followed the same notation used in the main text for the total number density, , the particle diameter, , and the molar composition, , where .

Apart from the HS contribution to the total free energy, we have an additional contribution due to the bonding free energy, , discussed in section II.C, which is a function of the probabilities for (see section II.C, Eqn. (7)). These probabilities are obtained through the law of mass action de las Heras et al. (2012) that in our case takes the form of two coupled non-linear equations:

(14) | |||||

(15) |

where all the interaction parameters needed for describing bonding between - and -patches enter in and de las Heras et al. (2012):

(16) |

(17) |

In our work we have approximated and by using the contact values of the partial radial distribution functions, , , and (where , for a binary mixture of hard spheres as obtained from the Percus-Yevick Equation Lebowitz and Rowlinson (1964):

(18) |

(19) |

where:

(20) |

The bonding volumes and present in Eqs. (16) and (17) are given by:

(21) |

(22) |

where and () are respectively the interaction ranges and the angular patch widths presented in sections II.A and II.B. With our choices for the interaction parameters we have . Once Eqs. (14) and (15) are solved (by using Eqs. (16)-(22)), the probabilities and that an -site is specifically bonded to another - or to a -site are obtained by the relations: and .

We finally present the procedure for constructing the phase diagram of our binary system, i.e locating its coexisting points, through the minimization of the Gibbs free energy once we have the analytical expression for the total Helmholtz free energy, . First we obtain, for a fixed composition , temperature , and for a given (target) pressure , those densities that satisfy the equation of state:

(23) |

In principle, we can have different densities which are solution of Eqn. (23) (all of them having a common pressure, composition, and temperature). From these densities we take that one which minimizes the Gibbs free energy per particle of our system, , i.e. we take that density for which our system is most stable. This process is then repeated for a new composition by keeping and constant, thus covering the whole range of compositions, . As a result we obtain the stable values of the Gibbs free energy for any composition (and therefore for any partial density, ) for a fixed pressure and temperature.

Figure 11 shows an example of at in which we have two points and that share a common tangent, that is:

(24) |

(25) |

In Eqs. (24) and (25) we recognize the conditions of chemical equilibrium of two phases with compositions and , corresponding to the total number densities and which are the solutions of Eqn. (23) for and . In other words, Eqs. (24) and (25) are equivalent to the equalities of the chemical potentials and of the two species and in the two coexisting phases given by the compositions and . That is:

(26) |

(27) |

where (see Fig. 11) :

(28) |

(29) |

with .

We should remember that Eqs. (24) and (25), as obtained from Eqn. (23), do not only consider the chemical equilibrium of the two species in the two phases given by and but also assure thermal and mechanical equilibrium since all the points with which we construct have a common temperature and pressure. Additionally, it is obvious that for those pressures and temperatures for which is convex ( ) we will have no phase coexistence. To construct the complete phase diagram corresponding to a fixed temperature (see Fig. 2) we repeat the whole process by varying .

Finally, tie lines connecting two coexisting points and (see Fig. 2) are defined by all those points that separate into two phases with compositions and and total number densities and conserving the total number of - and -particles as well as the total volume. These constraints result in a function (tie line) for any double pair and :

(30) |

with domain (here we assume without loss of generality) and satisfying:

(31) |

## References

- Wang et al. (2012) Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng, A. D. Hollingsworth, M. Weck, and D. J. Pine, Nature 491, 51 (2012).
- Jones and Mirkin (2012) M. R. Jones and C. A. Mirkin, Nature 491, 42 (2012).
- Chen et al. (2011) Q. Chen, S. C. Bae, and S. Granick, Nature 469, 381 (2011).
- Ruzicka et al. (2011) B. Ruzicka, E. Zaccarelli, L. Zulian, R. Angelini, M. Sztucki, A. MoussaÃ¯d, T. Narayanan, and F. Sciortino, Nature Mater. 10, 56 (2011).
- Seeman (2003) N. C. Seeman, Nature 421, 427 (2003).
- Condon (2006) A. Condon, Nature Rev. Genet. 7, 565 (2006).
- Stewart and McLaughlin (2004) K. M. Stewart and L. W. McLaughlin, J. Am. Chem. Soc. 126, 2050 (2004).
- Biffi et al. (2013) S. Biffi, R. Cerbino, F. Bomboi, E. M. Paraboschi, R. Asselta, F. Sciortino, and T. Bellini, Proceedings of the National Academy of Sciences 110 15633 (2013).
- Kraft et al. (2012) D. J. Kraft, R. Ni, F. Smallenburg, M. Hermes, K. Yoon, D. A. Weitz, A. van Blaaderen, J. Groenewold, M. Dijkstra, and W. K. Kegel, Proceedings of the National Academy of Sciences 109, 10787 (2012).
- Jiang and Granick (2012) S. Jiang and S. Granick, eds., Janus particle synthesis, self-assembly and applications, RSC Smart Materials (The Royal Society of Chemistry, 2012).
- Walther and Muller (2008) A. Walther and A. H. E. Muller, Soft Matter 4, 663 (2008).
- Chen et al. (2012) Q. Chen, S. C. Bae, and S. Granick, J. Am. Chem. Soc. 134, 11080 (2012).
- Y. Iwashita (2013) Y. K. Y. Iwashita, submitted (2013).
- Glotzer and Solomon (2007) S. C. Glotzer and M. J. Solomon, Nature Mater. 6, 557 (2007).
- de Graaf et al. (2011) J. de Graaf, R. van Roij, and M. Dijkstra, Phys. Rev. Lett. 107, 155501 (2011).
- Romano and Sciortino (2011) F. Romano and F. Sciortino, Nature Mater. 10, 171 (2011).
- Angioletti-Uberti et al. (2012) S. Angioletti-Uberti, B. M. Mognetti, and D. Frenkel, Nature Mater. 11, 518 (2012).
- Smallenburg et al. (2012) F. Smallenburg, L. Filion, M. Marechal, and M. Dijkstra, Proceedings of the National Academy of Sciences 109, 17886 (2012).
- Bianchi et al. (2011) E. Bianchi, R. Blaak, and C. N. Likos, Phys. Chem. Chem. Phys. 13, 6397 (2011).
- Preisler et al. (2013) Z. Preisler, T. Vissers, F. Smallenburg, G. Munao, and F. Sciortino, Journal of Phys. Chem. B 117, 9540 (2013).
- Hong et al. (2006) L. Hong, A. Cacciuto, E. Luijten, and S. Granick, Nano letters 6, 2510 (2006).
- Capone et al. (2012) B. Capone, I. Coluzza, F. LoVerso, C. N. Likos, and R. Blaak, Phys. Rev. Lett. 109, 238301 (2012).
- Zaccarelli (2007) E. Zaccarelli, J. Phys. Cond. Mat 19, 323101 (2007).
- Sciortino and Zaccarelli (2011) F. Sciortino and E. Zaccarelli, Current Opinion in Solid State and Materials Science 15, 246 (2011), ISSN 1359-0286.
- Sear and Gelbart (1999) R. P. Sear and W. M. Gelbart, J. Chem. Phys. 110, 4582 (1999).
- Groenewold and Kegel (2001) J. Groenewold and W. K. Kegel, J. Phys. Chem. B 105, 11702 (2001).
- Sciortino et al. (2004) F. Sciortino, S. Mossa, E. Zaccarelli, and P. Tartaglia, Phys. Rev. Lett. 93, 055701 (2004).
- Tlusty and Safran (2000) T. Tlusty and S. A. Safran, Science 290, 1328 (2000).
- Russo et al. (2011) J. Russo, J. M. Tavares, P. I. C. Teixeira, M. M. Telo da Gama, and F. Sciortino, Phys. Rev. Lett. 106, 085703 (2011).
- Roldan-Vargas et al. (2013) S. Roldan-Vargas, F. Smallenburg, W. Kob, and F. Sciortino, Scientific Reports 3, 2451 (2013).
- Pawar and Kretzschmar (2008) A. B. Pawar and I. Kretzschmar, Langmuir 24, 355 (2008).
- Romano and Sciortino (2012) F. Romano and F. Sciortino, Nature Commun. 3, 975 (2012).
- Sciortino et al. (2009) F. Sciortino, A. Giacometti, and G. Pastore, Phys. Rev. Lett. 103, 237801 (2009).
- Sciortino et al. (2010) F. Sciortino, A. Giacometti, and G. Pastore, Physical Chemistry Chemical Physics 12, 11869 (2010).
- Reinhardt et al. (2011) A. Reinhardt, A. J. Williamson, J. P. Doye, J. Carrete, L. M. Varela, and A. A. Louis, J. Chem. Phys. 134, 104905 (2011).
- Rovigatti et al. (2013a) L. Rovigatti, J. M. Tavares, and F. Sciortino, Phys. Rev. Lett. 111, 168302 (2013).
- Almarza (2012) N. G. Almarza, Phys. Rev. E 86, 030101 (2012).
- Binder and Kob (2011) K. Binder and W. Kob, Glassy materials and disordered solids (World Scientific, Singapore, 2011).
- Saika-Voivod et al. (2013) I. Saika-Voivod, F. Smallenburg, and F. Sciortino, Accepted in J. Chem. Phys., arXiv:1309.2198 (2013).
- Wertheim (1984a) M. S. Wertheim, J. Stat. Phys. 35, 19 (1984a).
- Wertheim (1984b) M. S. Wertheim, J. Stat. Phys. 35, 35 (1984b).
- Rovigatti et al. (2013b) L. Rovigatti, D. de las Heras, J. M. Tavares, M. M. T. da Gama, and F. Sciortino, J. Chem. Phys. 138, 164904 (2013b).
- Kern and Frenkel (2003) N. Kern and D. Frenkel, J. Chem. Phys. 118, 9882 (2003).
- Barcenas et al. (2007) M. Barcenas, J. Douda, and Y. Duda, J. Chem. Phys. 127, 114706 (2007).
- Rapaport (2009) D. C. Rapaport, Progress of Theoretical Physics Supplement 178, 5 (2009).
- de la Peña et al. (2007) L. H. de la Peña, R. van Zon, J. Schofield, and S. B. Opps, J. Chem. Phys 126, 074105 (2007).
- Smallenburg and Sciortino (2013) F. Smallenburg and F. Sciortino, Nature Phys. 9, 554 (2013).
- Virnau and Müller (2004) P. Virnau and M. Müller, J. Chem. Phys. 120, 10925 (2004).
- Mansoori et al. (1971) G. A. Mansoori, N. F. Carnahan, K. E. Starling, and J. T. W. Leland, J. Chem. Phys. 54, 1523 (1971).
- de las Heras et al. (2012) D. de las Heras, J. M. Tavares, and M. M. Telo da Gama, Soft Matter 8, 1785 (2012).
- Rovigatti and Sciortino (2011) L. Rovigatti and F. Sciortino, Mol. Phys. 109, 2889 (2011).
- de las Heras et al. (2011) D. de las Heras, J. M. Tavares, and M. M. T. da Gama, J. Chem. Phys. 134, 104904 (2011).
- Smallenburg et al. (2013) F. Smallenburg, L. Leibler, and F. Sciortino, Phys. Rev. Lett. 111, 188002 (2013).
- Montarnal et al. (2011) D. Montarnal, M. Capelot, F. Tournilhac, and L. Leibler, Science 334, 965 (2011).
- Foffi and Sciortino (2007) G. Foffi and F. Sciortino, J. Phys. Chem. B 111, 9702 (2007).
- Bianchi et al. (2006a) E. Bianchi, J. Largo, P. Tartaglia, E. Zaccarelli, and F. Sciortino, Phys. Rev. Lett. 97, 168301 (2006a).
- Flory (1953) P. J. Flory, Principles of Polymer Chemistry (Cornell University Press, 1953).
- Flory (1941) P. J. Flory, J. Am. Chem. Soc. 63, 3096 (1941).
- Romano et al. (2010) F. Romano, E. Sanz, and F. Sciortino, J. Chem. Phys. 132, 184501 (pages 9) (2010).
- Romano et al. (2011) F. Romano, E. Sanz, and F. Sciortino, J. Chem. Phys. 134, 174502 (2011).
- Kumar and Panagiotopoulos (2013) S. K. Kumar and A. Z. Panagiotopoulos, Private communication (2013).
- Foffi et al. (2002) G. Foffi, K. A. Dawson, S. V. Buldyrev, F. Sciortino, E. Zaccarelli, and P. Tartaglia, Phys. Rev. E 65, 050802 (2002).
- Zaccarelli et al. (2002) E. Zaccarelli, G. Foffi, K. A. Dawson, S. V. Buldyrev, F. Sciortino, and P. Tartaglia, Phys. Rev. E 66, 041402 (2002).
- Bianchi et al. (2008) E. Bianchi, P. Tartaglia, E. Zaccarelli, and F. Sciortino, J. Chem. Phys. 128, 144504 (2008).
- Lebowitz and Rowlinson (1964) J. L. Lebowitz and J. S. Rowlinson, J. Chem. Phys. 41, 133 (1964).