Gelling by heating.

We exploit the concept of competing interactions to design a binary mixture of patchy particles that forms a reversible gel upon heating. Our molecular dynamics computer simulation of such a system shows that with increasing temperature the relaxation dynamics slows down by more than four orders of magnitude and then speeds up again. The system is thus a fluid both at high and at low temperatures and a solid-like disordered open network structure at intermediate temperature. We further discuss the feasibility of realizing a real material with this reversible behavior.

Some of the most versatile and efficient strategies for designing new materials with unconventional behavior are based on the idea of competitive interaction. While in biology the term competitive interaction indicates the result of rivalry between two or more species competing for resources, in physics it stands for the presence of several interaction mechanisms that can stabilize competing local structures, leading to novel and highly interesting features of the system. Recent examples for this mechanism in soft matter systems include, among others, the competition between short-range attraction and long-range repulsion in charged colloids giving rise to cluster phases123, the competition between chaining and branching in patchy colloids45 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 gas5, and the design of DNA-coated colloids with two different DNA sequences for the purpose of establishing a competition between intra and inter-particle interactions, favoring crystal formation6. Often the very nature of these competing mechanisms promotes the emergence of a structure controlled by energy (stable at low temperature, T) which competes with a structure stabilized by entropy at intermediate T.

Recent progress in the synthesis of colloids has led to a novel type of particles that have highly directional and selective interactions, thus providing valence to colloids7. In particular the recent developments in large scale synthesis of patchy colloidal particles8 and in DNA nanotechnology have already allowed for the production of versatile constructs in which specificity, valence and interaction strength can be tuned at will910. This ability to control the interactions between nano-and meso-sized particles8, and to design the geometric properties of the patches1112 and/or their functionalization, offers now the possibility to take advantage of the competitive interactions that are usually present in such systems and hence to modulate material properties with external control parameters. In this Report we develop one such possibility: The design of a material that continuously (i.e. without the interference of a first order transition) and reversibly gels upon heating. In contrast to gelation on heating observed in protein solutions13, where kinetic arrest is driven by the irreversible denaturation of the proteins, here we require thermo-reversibility of the process and a simple design of the constituent particles. We show how a binary mixture of limited valence particles can provide a model in which the competition between entropy and potential energy causes the system to show a re-entrant behavior, passing reversibly from a fluid to a gel and again to a fluid if T is varied. The proposed model, amenable to theoretical analysis and experimental realization, shows in simulations a slowing down of the dynamics of several orders of magnitude on heating, providing the first example of controlled reversible gelation.

Patchy colloidal particles with a limited number f of attractive patches progressively cluster when cooled, leading such systems to form a percolating network that at sufficiently low T will contain all particles in the system141516171819. Due to this aggregation, the dynamics slows down by several orders of magnitude. In the network state the lifetime of the bonds (a T-controlled quantity) fixes the timescale over which the system behaves as a solid. It is known that if f = 4, the particles (here called A species) form a random tetrahedral network which closely resembles the structure of network-forming atomic systems like silica or silicon20. To melt the network at low T, we add a second species (B) with only a single patch which competes for bonding with the patches on the network forming A-species21. The idea is to design a competitive mechanism such that the bonding between A- and B-particles becomes dominant, but only at a temperature which is much lower than the one at which the AA-network is formed. As a result, the stable low-T phase consists of A-particles decorated with f B-particles which are free to diffuse in the sample volume, whereas at intermediate T the system forms a highly viscous AA-network that is progressively fragmented and transformed into a fluid upon heating.

Results

Model

We consider a binary mixture of patchy colloids where each A-particle has f = 4 patches on its surface that are arranged in a tetrahedral geometry and the B-particles have only one patch (see Fig. 1). The patch-patch interaction is modeled via a Kern-Frenkel potential22, a model that has been extensively used over the last decade to compare simulations and experiments on the self-assembly of patchy colloids232425. Each A-patch can interact either with a B-patch with unit energy and bonding volume or with another A-patch with energy and . No BB-bonding is allowed. The bonding volumes are each determined by an interaction range δ and an angular patch width θ (α ∈ {AA,AB}) (see Fig. 1 and Methods). The attractive patch-patch interaction is complemented by an isotropic hard-core repulsion, where the spherical cores have diameters σ and σ = 0.35σ. The size ratio was chosen such that the B-particles can block the A-patches from bonding to other A-patches without significantly increasing the packing fraction of the pure A system. Due to geometric constraints, each patch can be involved in only a single bond.

Figure 1

The Kern-Frenkel model.

(a) Schematic of the interaction parameters in the Kern-Frenkel model. An A-patch can bond with either another A-particle (b), or with a B-particle (c). Panel (d) shows a “flower”, i.e. a fully bonded cluster consisting of one A-particle and four B-particles, representing the lowest-energy state of the system. Here the interaction ranges and the angular patch widths are cos θ = 0.92, δ = 0.15σ, and cos θ = 0.99, δ = 0.2σ. With these choices the bonding volumes are given by and (see Methods). Note that in the binary Kern-Frenkel model, θ and δ of each bond are defined by the species of both bonding partners and are not properties of individual particles.

Using event-driven molecular dynamics simulations262728 (see Methods), we have studied a system of N = 600 and N = 2400 particles, corresponding to a total number density , with partial number densities and , for a wide range of T, whose unit is given by , where k is Boltzmann's constant. The composition of the system is thus fixed at x = 0.2. The density of A-particles corresponds to the optimal density at which tetrahedral particles form an unstrained fully bonded network29. With this composition, the fully bonded network has an energy of , whereas a configuration in which all the B-particles are bonded to the A-particles has the significantly lower energy of . We also simulate for a low a reference system composed of 600 flowers (see Fig. 1d), i.e. A-particles bonded to four B-particles. At this low T no bond breaking events take place within the simulation time.

Temperature dependence of the bonding probability

Figure 2 demonstrates the effect of the competitive interactions present in our system. It shows the T–dependence of the probability that a patch on an A-particle is bonded to another A-patch (pAA) or to a B-patch (pAB). On cooling, pAA starts to grow, signaling the onset of the network formation, and reaches a maximum around . We recall that within a mean-field description, percolation of particles with valence f = 4 takes place at pAA = 1/330. Since at the maximum we find pAA ≈ 0.9 we can conclude that at this T the A-particles have formed a highly bonded percolating network. Upon further cooling, pAA starts to decrease whereas pAB grows quickly, showing that the AB-bonds are starting to replace the AA-bonds. The entropy associated with the larger bonding volume for the AA interaction is crucial for promoting the formation of a large number of AA-bonds at intermediate T, before the energetically preferred but entropically disfavored AB-bonds set in. Figure 2 also shows the parameter-free theoretical predictions for the bonding probabilities as obtained from the first-order thermodynamic perturbation theory developed by Wertheim313233 (details on the Wertheim calculations are reported in the Methods section). The Wertheim theory nicely captures the mechanism of competing interactions, reproducing the position and height of the maximum of pAA as well as the low T trends of pAA and pAB.

Figure 2

Bonding probabilities.

Probability that a patch on an A-particle is bonded to another patch on an A-particle, pAA (black circles), or to a patch on a B-particle, pAB (red squares). The dashed lines are the prediction for these probabilities as calculated from the Wertheim theory.

Structure

Figure 3 shows the unusual T dependence of the structure of the system, related to the non-monotonic behavior of p. At high T, the partial structure factor S(q) shows the conventional q–dependence found in simple liquids with a main peak around qσ = 7.2. Upon decreasing T the main peak splits into two, one located around qσ = 5.2 and a higher one around 8.4. This double peak feature is typical of liquids that have a local tetrahedral network structure, such as silicon or silica20. The peak at qσ ≈ 8.4 corresponds to the nearest neighbor distance between two bonded A-particles, whereas the one at around 5.2 is associated with the second-nearest neighbors in the tetrahedral network. Note that this double peak structure is most pronounced at , i.e. at the T at which p has a maximum (see Fig. 2) and hence the gel is maximally connected. When T is lowered even further the double peak structure disappears and S(q) becomes again similar to the structure factor of a fluid composed of flowers (which is also included in Fig. 3 as a reference). Since the size of a flower is larger than that of an A-particle, the peak position at low T is to the left of the one observed at high T.

Figure 3

Structure factor.

Partial structure factor SAA(q) for different values of T (solid lines with different symbols). The SAA(q) for a fluid of flowers at (black solid line) is also represented.

Temperature dependence of the diffusion coefficient

We now quantify the effect of the competing interactions on the dynamics of the system and provide evidence that the change of the structure rich in AB-bonds at low T to the highly bonded AA-network generates a slowing down of the dynamics on heating. For this, we calculate the mean squared displacement (MSD) for the particles of both species and then their corresponding diffusion coefficients D (α ∈ {A,B}) from the long-time behavior of the MSD via the Einstein relation. We should note that the center of mass (CM) of a single species has a non-zero velocity (which is compensated by the CM motion of the other species). To obtain meaningful results for the MSD of a single species, we subtract the CM drift of the species in question before evaluating the MSD. In addition, to remove the trivial trend originated from the T-dependence of the thermal velocity we divide D by a reference diffusion coefficient , where and m is the mass of an A-particle. Figure 4 shows the T-dependence of D in an Arrhenius plot. At high T, D is approximately constant for both type of particles, indicating that bonds do not play a significant role. On cooling, D starts to decrease very rapidly, with a super-Arrhenius T-dependence reminiscent of that observed in molecular networks20, turning into an Arrhenius law with an activation energy approximately equal to (see dashed line in Fig. 4). Similar values of the activation energy are typically found in tetrahedral network-forming systems where most of the particles belong to the percolating cluster, and bond breaking is the bottleneck for relaxation183435. Before the gel starts to decompose at a temperature below , D has already decreased by four orders of magnitude compared to its value at high T, indicating the formation of a persistent network. For , D starts to increase. This rising persists down to the lowest T at which we were able to equilibrate the system. We emphasize that this non-monotonic T-dependence is only observed for the A-particles, i.e. the particles which are involved in the formation of the network. In contrast, D shows only a rather mild T-dependence. Fig. 4 also shows the diffusion coefficient of the fluid of flowers at for which D = D, the common value to which both D and D converge.

Figure 4

Diffusion coefficient.

Arrhenius plot of the normalized diffusion coefficient D/D0 (α = A, circles, α = B, squares). Also included is an Arrhenius law with activation energy (dashed line). The diffusion coefficient of a fluid of flowers at low T is represented by a blue diamond.

Relaxation times

A non-monotonic behavior of the characteristic time is also found in the time evolution of the collective- and self-intermediate scattering functions. Their study provides insight into how the relaxation dynamics depends on the considered length scale. Figure 5a shows an Arrhenius plot of the relaxation time τ(q) determined from the time integral of the intermediate scattering function of the A-particles for two different q–vectors: qσ = 5.2 and qσ = 7.2, which correspond, respectively, to the location of the first peak in the network and in the high-T fluid (see Fig. 3). We find that the self and collective relaxation times τ(q), normalized by τ0, show qualitatively the same T-dependence: a plateau at high T, a fast increase within the T-range in which the network is formed, a quick decrease once the network starts to break up again, and a final plateau at low T. This T-dependence is observed for both values of q, indicating that the relaxation mechanism does not depend on the length scale considered.

Figure 5

Relaxation times.

(a) Arrhenius plot of the normalized relaxation time τ(q)/τ0 as obtained from the self (dashed lines) and collective (solid lines) scattering functions, where . The different curves correspond to the wave-vectors given by the first two peaks in S(q). (b) Arrhenius plot of the normalized bond-persistence time τ/τ0 for the AA (squares) and AB bonds (circles). Also included is an Arrhenius law with activation energy (blue solid line).

To provide further evidence that the system is ergodic on long time scales, i.e. that the structure of the system has completely lost its memory of the initial state, we have investigated the bond persistence function p(t), i.e. the probability that a bond which is present at time zero is also present at time t (see Supplementary Information). Hence, the relaxation time of p(t) provides information on the restructuring time of the network connectivity. Figure 5b shows an Arrhenius plot with the T-dependence of the decay time τ, where p(τ) = e−1. At intermediate and low T, τ is larger than the relaxation times shown in Fig. 5a. We find thus that p(t) decays to zero only on a time scale that is significantly longer than the relaxation times associated with the scattering functions, confirming that some fraction of spacial decorrelation of the network, as quantified by the collective scattering function, takes place at partially fixed bonding pattern. The extreme case occurs at very low T where the system is a fluid of flowers that relaxes relatively quickly but in which AB-bonds survive for a very long time. Figure 5b also shows that the T-dependence of the lifetime of the AA-bonds differs strongly from the one of the AB-bonds. The latter shows in the whole T–range an Arrhenius dependence with an activation energy very close to (see blue solid line in Fig. 5b). This suggests that the mechanism for the breaking of an AB-bond is not collective in nature but a simple activated process. In contrast, the breaking time for an AA-bond follows an Arrhenius law at high T but becomes super-Arrhenius within the T-region in which the gel forms, showing that this bond-breaking process is of collective nature. At very low T the effective activation energy becomes again an Arrhenius behavior with an activation energy given by .

Discussion

We have demonstrated that a judicious choice of the interaction parameters of a binary mixture of A and B patchy particles shows a non-monotonic and reversible T-dependence of its dynamic properties. We set up a competition between network-forming AA-bonds and network-breaking AB-bonds, which are favored by entropy and energy, respectively. Thus, at intermediate T, entropy stabilizes a viscous network gel, whereas both at low and high T, the network breaks up, leading to a simple fluid state. Our results are robust with respect to a change in the model parameters and we have indeed observed a qualitatively similar phenomenology in our system for a wide range of the interaction energy and density, and . Since our system requires no fine-tuning, we expect that the observed phenomenon can also be realized experimentally. In particular, two experimental systems are very promising candidates to test the ideas presented in this Report: (i) a solution of DNA constructs of valence four3637 in the presence of competing DNA single strands and (ii) the binary mixture of patchy particles recently synthesized in Ref. 8. Both of these systems offer single-patch and four-patch colloidal particles with controllable interaction strengths that are already amenable for experimentation, and thus have the potential to provide soft materials that reversibly gel upon heating.

Methods

Kern-Frenkel model

The particles investigated in the paper are modeled via the well-known Kern-Frenkel model22. In this model, particles interact by means of a combination of a hard-sphere potential uHS and an attractive directional interaction upatch. The hard-core repulsion between two particles i and j is given by: here, β = 1/k, with k Boltzmann's constant and T the temperature, r is the center-to-center distance between the particles, and σ = (σ + σ)/2 denotes the contact distance between the particles, with σ(σ) the diameter of particle i(j). 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 when 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 (see Fig. 1 of main text). The size of the patches is determined by an opening angle θ. The potential energy of two particles i and j is thus given by where uSW is a square-well attraction, given by: The function Φ(r, {p}) is defined as where r = r − r, {p} is a set of normalized vectors pointing from the center of particle i towards the center of each of its patches, and is a unit vector in the direction of r.

In the model used in the Report, all particles are either species A or species B. The A-particles (with hard sphere diameter σ) have four patches each, arranged in a tetrahedral geometry. The B-particles are smaller, with σ = 0.35σ, and have only a single patch. The parameters for the interaction between two A-particles are given by cos θ = 0.92, and δ = 0.15σ. For the AB-interactions, cos θ = 0.99, and δ = 0.2σ. There are no attractive interactions between the B-particles. The ratio between the two interaction strengths is fixed at .

Event-driven molecular dynamics

To study the relaxation dynamics of our model we use event-driven Molecular Dynamics (EDMD) simulations2627. The implementation of the EDMD simulation relies on the numerical prediction of bond formation and bond breaking events, and follows the same scheme as described in Ref. 28. In the model under consideration here, the mass m 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 were also 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 velocity and angular velocity, drawn from a Maxwell-Boltzmann distribution. While measuring the diffusion coefficient, no thermostat is used, so that the total energy in the system is constant. To speed up equilibration of the system at low temperature, the initial configuration is taken from a standard canonical Monte Carlo simulation equilibrated at the same temperature T.

Wertheim theory

Wertheim's thermodynamic perturbation theory allows us to obtain an analytical expression for the free energy of pure fluids and fluid mixtures (a detailed description can be found in Refs. 31, 32). In the context of Wertheim's theory, and specialized to our system, the probability p that an α-site (α ∈ {A, B}) is bonded is obtained through the law of mass action, which in our case takes the form of a set of two coupled equations33: where x = 0.2 is the molar fraction of the species A, ρ is the total number density. All interaction parameters needed for describing bonding between AA-and AB enter in Δ and Δ: where and are the bonding interaction energies of the AA- and AB-bonds. In the present work we have approximated Δ and Δ by using the contact values of the partial radial distribution functions33, g(σ), g(σ), and g(σ) (where σ = (σ + σ)/2), for a binary mixture of hard spheres as obtained from the Percus-Yevick Equation38: where , , being ρ and ρ the partial number densities of the different species. The bonding volumes and present in Equations (7) and (8) are given by: where δ and θ (γ ∈ {AA, AB}) are respectively the interaction ranges and the angular patch widths defined in the model section. Once Eqs. (5) and (6) are solved (by using Eqs. (7)–(12)), the probabilities pAA and pAB that an A-site is specifically bonded to another A- or to a B-site are obtained by the relations: pAA = p − p and pAB = p.

Author Contributions

F. Smallenburg and F. Sciortino wrote the code for the EDMD and MC simulations. S. R.-V. and W. K. performed the simulations and the analysis of the data. All the authors contributed to the interpretation of the results and the writing of the paper.