Main

The ability to control the selectivity and angular flexibility of interparticle interactions5 makes it possible to provide valence to colloids. This can be done by means of chemical6,7,8 or physical9 patterning of the colloid surface, that is, locally changing either its chemical properties (for example, hydrophobicity), or the physical properties (for example, roughness). Exploiting valence offers enormous possibilities, some of which have been addressed in recent years, theoretically10, numerically4,11,12 and experimentally2,7,13,14. As colloidal interactions are typically short-ranged, valence provides an implicit quantization of the particle energy, which becomes essentially proportional to the (limited) number of formed bonds. In an optimal arrangement, all particles bind with their maximum number of neighbours and the system is in its lowest possible (ground) energy state. Typically, such an optimal arrangement is spatially ordered, defining the most stable crystal phase(s).

In principle, when the density is not too high, it is possible to imagine disordered arrangements in which all particles are fully bonded, with exactly the same energy as the crystal. Such a state, which was recently shown to be stable in a model system for DNA-functionalized colloids12, is a liquid: a fluid, disordered phase at temperatures below the critical temperature. Under this unconventional condition, the stability of the system becomes controlled by its entropy, a case reminiscent of (purely entropy-driven) hard-sphere particles. In this study we demonstrate that the flexibility of the bonds (encoded in the angular patch width)—a tunable quantity in the design of patchy colloidal particles—is the key element in controlling the entropy of the liquid as compared with that of the crystal. Large binding angles combined with limited valence give rise to a thermodynamically stable, fully bonded liquid phase.

The model we consider is an extension of the widely used Kern–Frenkel model15 for patchy particles, where we explicitly enforce a single-bond-per-patch condition. Each spherical particle of diameter σ has a number f of circular attractive patches on its surface, characterized by a maximum opening angle θmax. Two particles are bonded with a fixed bonding energy ε when their centre-to-centre distance is smaller than the interaction range σ+δ, and the vector connecting the particles passes through a patch on both particles (see Supplementary Information). We examine monodisperse systems of particles with f = 4 attractive patches, arranged in a tetrahedral configuration. In addition, to show that the behaviour of these systems is robust with respect to changes in the patch positions, as well as the number of patches, we also examine (for a single patch size) the phase diagram of particles with f = 3 patches, arranged along the equator. Although the one-bond-per-patch limitations renders this model inadequate for colloidal particles with large attractive surface regions, we stress here that the model mimics systems with a limited-valence bonding mechanism that retains a large degree of bonding flexibility for all experimentally relevant temperature scales. For example, one might consider colloidal particles functionalized with a limited number (f) of long strands of DNA (refs 12, 16), DNA constructs that can form a fixed number of bonds14,17, or even polymer networks18 where the particles are polymerizable monomers with a fixed functionality.

The phase diagrams, evaluated using free-energy calculations (see Methods), for different θmax are shown in Fig. 1 as a function of density ρ and T. The stable phases include the disordered gas and liquid phases, a diamond cubic crystal phase, a body-centred cubic (bcc) crystal phase, consisting of two interlocking diamond lattices, and a face-centred cubic (fcc) crystal phase, where, for sufficiently narrow patches, the bonds are periodically ordered at low T and disordered at high T. All crystals are fully bonded at low T.

Figure 1: Phase diagram of patchy colloids for different patch widths.
figure 1

Phase diagrams in the dimensionless density (ρ σ3)–temperature kBT/ε representation, where kB is the Boltzmann constant, for patchy particles. a–e, Diagrams for particles with four tetrahedrally arranged patches and patch widths cosθmax = 0.9 (a), cosθmax = 0.8 (b), cosθmax = 0.7 (c) and cosθmax = 0.6 (d), as well as a diagram for particles with three patches along the equator (e) and cosθmax = 0.7. The interaction range in all cases is fixed at δ = 0.12σ, a value typical of colloidal interactions. For tetrahedrally arranged patches, narrower patches (cosθmax>0.9) behave as the cosθmax = 0.9 case19 and are not shown here. The grey areas denote coexistence regions and the points denote calculated coexisting states (tie lines are horizontal). Black triangles indicate triple points. The points below kBT/ε = 0.1 are based on extrapolation of the potential energy. The label ofcc denotes the orientationally ordered fcc phase; the disordered fcc phase is simply denoted as fcc. f–i, Pictures showing the unit cell of bcc (f), with the two interspersed diamond lattices indicated in different colours, diamond (g), orientationally ordered fcc (h) and disordered fcc (i).

Full size image

For narrow patch width, the four-patch phase diagram closely resembles the standard phase diagrams with a triple point below which the liquid state ceases to exist19. Unexpectedly, for wider patches the diamond crystal disappears from the phase diagram and a region in which the liquid is stable down to vanishing T opens up at intermediate ρ, in agreement with a numerical study of DNA-functionalized colloids12,20. In this range of densities, the liquid is the thermodynamically stable phase, despite its intrinsic disorder. Increasing the patch width even further reduces the region of stability for the crystal phases. For very wide patches, the phase diagram simply consists of gas, liquid and fcc crystal phases, with a large stable liquid region in the zero-temperature limit. A typical snapshot of a fully bonded liquid is shown in the Supplementary Information. For the three-patch case, the phase diagram for narrow patches is known to contain an fcc crystal phase at high ρ, and a crystal of interpenetrating hexagonal sheets or separate hexagonal sheets at intermediate ρ (ref. 21). For wider patches (cosθmax = 0.7), we see again that the lower-density ordered phases disappear in favour of a fully bonded liquid phase, leaving disordered fcc as the only stable crystal structure. The remainder of this paper will focus on the particles with four patches.

Why is the liquid more stable than the crystal? Figure 2 shows the T dependence of the number of bonds per particle nb for different θmax. On cooling, nb increases continuously to four, and the system progressively approaches the fully bonded random tetrahedral network state, the ground state of the system. At low T, the potential energy of the liquid is thus equal to that of the crystals. Hence, the stability of the liquid, as in the hard-sphere case, results from a subtle competition between vibrational Svib and configurational Sconf components of the total entropy Stot = Svib+Sconf (refs 3, 22), both of which can be calculated (see Methods). Svib measures the (phase-space) volume explored by each particle in a fixed bonding topology, and is larger in the ordered structure than in the fluid. Sconf at T = 0 measures the number of distinct configurations resulting in a fully bonded macroscopic state. It is positive in the fluid phase, but vanishes in our system for the crystal phase. Note that Sconf can still play a role in crystals with asymmetric bonds, such as hydrogen bonds in ice23, or patchy particles with multiple patch types24.

Figure 2: Number of bonds per particle in the tetrahedral liquid phase.
figure 2

Here ρ σ3 = 0.57, the same density as the diamond crystal phase. In the zero-temperature limit, the number of bonds per particle nb reaches the maximum of f = 4. Points denote measurements in N V T simulations, and lines are fits to guide the eye. Note that the approach to the ground state (nb = 4) takes place at higher T for large θmax values, a result induced by the larger bonding volume.

Full size image

Figure 3a shows how patch width affects the entropy. As for the dense hard-sphere case, Svib in the diamond structure (SvibDC) is always larger than in the liquid (Svibliq). However, the entropy of the liquid is enhanced by Sconf, causing for wide angles the unconventional stability of the liquid phase even when T→0. In this large θmax region, the model provides a neat example of a system for which the Kauzmann temperature (TK; ref. 3) does not exist: Sconf (Fig. 3b) does not vanish when T→0. Thus, for wide patches, a previously unexplored thermodynamically stable state of matter arises: the disordered fully bonded network. Note that although the angular dependence of our bonding potential is flat, Sconf is large enough to stabilize the liquid even in the presence of a weak bending cost. This flexibility is realistic for single-strand DNA and polymer bonds, which have finite persistence lengths at all experimentally accessible temperature scales. In this case, the stability of the liquid phase may technically not hold all the way down to T = 0, but the stability of the liquid phase will still occur at any relevant T.

Figure 3: Total, vibrational and configurational entropy for the zero-temperature phases.
figure 3

Here ρ σ3 = 0.57, and each particle has four patches. a, The vibrational and total entropy for the diamond crystal (DC) and liquid (liq) phases. For the diamond crystal, the total entropy StotDCcoincides with the vibrational entropy SvibDC. Below cosθmax≈0.89, the fully bonded liquid becomes more stable than the crystal. For each angular patch width θmax, five independent Svibliq calculations are shown. Lines are polynomial fits to the points. With extremely lengthy simulations fully bonded configurations can be sampled up to cosθmax = 0.92. Note that the entropies are negative because we are investigating a classical system and the (constant) kinetic contribution is not included. b, Configurational entropy Sliqconf of the liquid at T→0, calculated as the difference between the fits for Stotliq and Svibliq.

Full size image

To provide further evidence that sampling of the fully bonded state is not pre-empted by dynamic arrest and to characterize the low- T dynamics we investigate the microscopic mechanism associated with the evolution of the network. We start by focusing on the defects of the fully bonded network (broken bonds), which act as elementary diffusing units. At any T close to zero the number of broken bonds per particle (αbb) depends on T as exp(−ε/2kBT) (see Supplementary Information). Figure 4a shows that αbb indeed satisfies the exponential T dependence for all θmax.

Figure 4: Defects in the tetrahedral network and their dynamic role.
figure 4

a, The number of broken bonds per particle αbb as a function of the inverse temperature ε/kBT, for different patch widths, at constant density ρ σ3 = 0.65 and f = 4. The line has slope −1/2. b, The (dimensionless) diffusion coefficient D τ/σ2 as a function of the inverse temperature ε/kBT. In the EDMD, an attempt to break, form or switch bond is performed with a rate γ. Data here refer to a constant bond switching rate γ τ = 100, a value for which the dynamics is not affected by the specific value of γ (see Supplementary Information). The lines have slopes −1/2 (for cosθmax = 0.6,0.7,0.8) and slope −2 (for cosθmax = 0.9).

Full size image

To quantify the low- T dynamics, we perform event-driven molecular dynamics (EDMD) simulations and examine the diffusion coefficient D at low T (Fig. 4b). We observe two different mechanisms for restructuring the network at low T: bond-breaking and bond-switching. For small θmax, network rearrangement requires the breaking and reforming of several bonds. As a result, D follows an Arrhenius law with an activation energy of approximately 2ε. However, for large θmax the bonding volumes of different neighbours can overlap and the switching of bonds with no energetic penalty allows the system to relax stresses much more quickly. In this case, D follows an Arrhenius law with an activation energy of only ε/2; the same as αbb. Thus, D is proportional to the number of broken bonds, demonstrating that these defects are the key element in the microscopic dynamics. The bond-switching mechanism mimics the microscopic dynamics of vitrimers18, a recently invented malleable network plastic where a catalyst enables the switching (transesterification) of bonds between nearby polymers. Interestingly, D always vanishes following the Arrhenius behaviour characteristic of strong atomic and molecular network-forming liquids25, regardless of the microscopic dynamics.

In summary, our calculations clearly show that for patchy colloids with a limited number of flexible bonds, the liquid smoothly forms a fully bonded disordered network on lowering T, through the progressive reduction of the number of network defects, without the intervention of a more stable crystal phase or phase separation. Counterintuitively, the low- T phase behaviour is dominated by entropy, rather than energy, as all competing phases reach the ground state. The main requirements for the existence of a stable liquid at low T are thus a large flexibility of the interparticle bonds (in this temperature range) and a low valence. The large patch width ensures that a wide variety of network realizations can be formed. This configurational entropy is instrumental in stabilizing the liquid phase with respect to the crystal. The low valence ensures that the density of the liquid coexisting with the gas is small10, leaving a large density region where the network can form. Both are key ingredients for the formation of what can be considered a previously unexplored state of matter: the thermodynamically stable fully bonded disordered network state.

Methods

Model.

We consider an extension of the Kern–Frenkel model15, where each patch can form at most one bond. In the normal Kern–Frenkel model, two patchy particles i and j, located at ri and rj respectively, feel an attraction given by

where ri j = rj−ri, {pi} is a set of normalized vectors pointing from the centre of particle i towards each of its patches, and usw is a square-well potential of hard diameter σ, range δ and depth ε.

Here, β = 1/kBT, with kB being Boltzmann’s constant.

The function Φ(rk l,{pk}) is defined as

In short, two particles bond if the vector connecting the centres of the particles passes through an attractive patch on each of the surfaces of each particle. Hence, as long as the patches on the same particle do not overlap, the Kern–Frenkel model allows for only a single bond between two particles. However, a single patch can make a bond with multiple different nearby particles if the patch width is larger (for our δ) than 0.895.

Here, we propose and implement (see Supplementary Information) a modification to the Kern–Frenkel model that allows only one bond per patch for all θmax values. In this modification, the overlap of the bonding volumes is not a sufficient criterion for defining the presence of a bond.

Simulation methods.

To calculate the free energy of the various phases as a function of the density, we use thermodynamic integration over the equation of state26, as measured using EDMD simulations. For the gas, the reference state is an ideal gas. For the liquid, and fluids above the critical point, the hard-sphere fluid is used as a reference state, integrating over the well depth ε at fixed density. Similarly, the hard-sphere crystal was used as a reference state for the disordered fcc crystal. We use the Frenkel–Ladd method to calculate free energies of the diamond, bcc and ordered fcc crystal structures27, using an additional aligning potential to fix the orientation of the particles28. Finally, we determine phase coexistence using a common-tangent construction, and use the resulting coexistence points to draw the phase diagram.

To determine which crystal structures to consider in our free-energy calculations, we employed the floppy-box method29, where candidate structures are generated by a Monte Carlo simulation with very few particles and a variable box shape at constant pressure. Resulting unit cells were then expanded to a larger system size and equilibrated in Monte Carlo N P T simulations. The structures that repeatedly formed and were still crystalline after this step were included in our free-energy calculations. For more details, see Supplementary Information.

For the calculation of the total entropy Stot at T = 0, we use thermodynamic integration over the temperature. As the potential energy of all phases decreases exponentially at low temperatures, extrapolation to zero temperature is straightforward, and integrating over the exponential decay allows us to directly calculate the entropy at T = 0. The common tangent constructions at zero temperature are shown in the Supplementary Information.

To calculate the vibrational entropy Svib in the liquid phase, we use the Frenkel–Ladd method on a fully bonded configuration taken from a simulation at low temperature. The network topology is kept fixed during the calculation. To reduce the maximum spring constant required in the integration, we first optimize the configuration in a separate MC simulation biased towards configurations where small displacements of particles do not cause overlaps or break bonds. See Supplementary Information for more information.

We use EDMD simulations for calculating the equations of state and potential energies required for thermodynamic integration, and to investigate the dynamics of the low-temperature liquid phase30,31. For a full description of the implementation of the Kern–Frenkel model in EDMD simulations, see Supplementary Information. For all EDMD simulations used here, the simulation box is chosen to be cubic or rectangular, with periodic boundary conditions. The simulations are performed at fixed number of particles N, volume V and temperature T. The temperature is regulated by means of a thermostat: at regular intervals, a random particle is given a new velocity and angular velocity, drawn from a Maxwell–Boltzmann distribution. Time is measured in units of , with m being the mass of a particle. For the simulations in this study, we chose the moment of inertia I = m σ2. Note that the choice of mass or moment of inertia has no effect on the equilibrium phase behaviour. Bond switching moves in the EDMD simulation (see Supplementary Information) are events that happen at a fixed rate γ for each patch in the system, such that per time unit τ each patch experiences γ τ attempts to form, switch or break a bond. Diffusion coefficients were obtained from the slope of the mean squared displacement as a function of time. The system sizes (ranging from N = 512 to 64,000) were chosen to ensure that the number of broken bonds in the system would be at least one on average, for the temperatures under consideration.