Directed percolation identified as equilibrium pre-transition towards non-equilibrium arrested gel states

The macroscopic properties of gels arise from their slow dynamics and load-bearing network structure, which are exploited by nature and in numerous industrial products. However, a link between these structural and dynamical properties has remained elusive. Here we present confocal microscopy experiments and simulations of gel-forming colloid–polymer mixtures. They reveal that gel formation is preceded by continuous and directed percolation. Both transitions lead to system-spanning networks, but only directed percolation results in extremely slow dynamics, ageing and a shrinking of the gel that resembles synaeresis. Therefore, dynamical arrest in gels is found to be linked to a structural transition, namely directed percolation, which is quantitatively associated with the mean number of bonded neighbours. Directed percolation denotes a universality class of transitions. Our study hence connects gel formation to a well-developed theoretical framework, which now can be exploited to achieve a detailed understanding of arrested gels.

Model systems of physical gels have been established, for example, dispersions of spherical colloidal particles, in which attractive forces are induced by a temperature variation or depletion effects due to the addition of linear non-adsorbing polymer chains 30,31 . Nevertheless, there are only few studies on the relation between the slowdown of the dynamics and the microscopic structure.
Here we show that the slowdown of the particle dynamics and the onset of ageing can be related to the microscopic structure of colloidal gels, namely a directed percolation (DP) transition. This is demonstrated in experiments on charged colloid-polymer mixtures 30,33 and in Brownian dynamics simulations of particles interacting through a combination of coulombic repulsion and depletion-induced attraction modelled by the Asakura-Oosawa potential 34 . The onset of DP is investigated along different paths in which either the attractive or the repulsive parts of the interactions are varied. We find that in any case continuous percolation precedes DP and that both can be associated with characteristic values of the mean number of bonds per particle. Furthermore, DP is found to be linked to the onset of ageing effects.

Results
States of the system. In charged colloid-polymer mixtures, two particle-particle interactions compete on different length scales: First, screened electrostatic interactions, which are long ranged for low salt concentrations and become shorter ranged on addition of salt; second, depletion interactions, which are short ranged, typically at most one tenth of the colloidal particle diameter, and controlled through the radius of gyration of the polymers. We experimentally investigated charged colloidal polymethylmethacrylate (PMMA) spheres with volume fraction FE0.2 and diameter s ¼ 1.72 mm in the presence of nonadsorbing linear polystyrene with radius of gyration r g ¼ 65 nm, where the effective polymer-colloid size ratio x eff E0.03 (Supplementary Table 1).
In experiments and simulations, we consistently observe four states that differ in the structural arrangements of the colloidal particles ( Fig. 1, left). Which state occurs depends on the salt concentration c salt , and therefore the screening length k À 1 characterizing the range of the repulsion and the polymer concentration c p , that is, attraction strength. We consider fixed c p and vary c salt , which will be referred to as 'path B'. In the absence of salt, mainly individual particles but also very few clusters are distributed homogeneously throughout the sample (sample C1; note that c p of this sample is higher than the one of the other samples, Supplementary Table 1). Adding salt, particles aggregate into isolated small clusters (sample B1). At large salt concentrations, larger clusters are observed which form a network in three dimensions (sample B2). Furthermore, a heterogeneous network structure with thick strands consisting of even larger, dense clusters are formed at the largest salt concentration studied (sample B3).
These states show different pair correlation functions g(r) (Fig. 1, right). Without salt (C1), g(r) is dominated by a peak located at the mean particle distance r À 1/3 . Furthermore, a small peak at particle contact reflects the presence of a small number of aggregated particles. For small amounts of added salt (B1), the increase in the contact peak and a corresponding decrease of the second peak of g(r) is consistent with the presence of doublets, triplets or small string-like clusters (visible in Fig. 1a). The size of the clusters is limited by the repulsive contribution to the potential 22,35,36 . For sample B2, the large peak at contact reflects the large fraction of particles forming clusters. Its g(r) is reminiscent of a percolated gel-like network structure 32,37,38 . Finally, the radial distribution function of sample B3 with its very high first peak at contact reflects the presence of a large number of bonded particles. Furthermore, its deep first minimum and pronounced split second peak can arise from triplet structures and local close packed arrangements. These features of the g(r) are observed in the experiments as well as in the simulations. While there is qualitative agreement between experimental and simulation results, small quantitative differences indicate that the potential and its parameter values used in the simulations are not perfectly describing the interactions present in the samples. (The mapping of the simulation and experimental parameters is described in Supplementary Note 1.) Nevertheless, both, experiments and simulations, reveal different states on the path towards a gel; homogeneously distributed individual particles, small clusters, networks with thin and thicker strands consisting of larger, more compact and denser clusters, respectively.
Characterization of the different states by number of bonds. The different structures are distinguished by the degree of particle aggregation. This can be linked to the mean number of bonds per particle, hNi. Individual particle-particle bonds cannot be determined unambiguously because the first peak in g(r) is broadened due to polydispersity (about 7%). Hence, two particles are defined to be bonded if they are closer than the position of the first minimum in g(r), which is located between 1.10 and 1.22s (Fig. 1, right).
The mean number of bonds per particle, hNi, is determined as a function of the interactions, which depend on the salt concentration c salt , determining the screening length k À 1 , and the polymer concentration c p , controlling the depletion part of the potential V D,min ¼ V D (s eff ) evaluated at the overall potential minimum V min ¼ V(s eff ) (Fig. 2). At low polymer, c p , and/or salt, c salt , concentrations, only very few bonds are found (white to purple), whereas high c p and/or c salt result in a bonded state (red). In contrast to previous results 23 , hNi monotonically grows without a re-entrant transition. Even for large screening and strong attraction, hNi hardly exceeds seven and hence the highly bonded state appears not to be in equilibrium. For a purely repulsive interaction 39,40 or in spring lattice models 41 such a change is expected once the system becomes isostatic, that is, globally stable with six bonds per particle.
The structures along two paths are analysed in more detail (Fig. 3). Along path A, the salt concentration c salt and hence the screening length k À 1 is kept constant while the polymer concentration c p is increased and hence the value V D,min of the potential at its minimum is decreased, that is, the attraction strength is increased. (Note that path A implies a high salt concentration, which cannot be achieved in the organic solvent mixture used in the experiments.) Instead, along path B, the screening length k À 1 is reduced at constant attraction strength. On increasing attraction on path A, hNi increases and the distribution of the number of bonds per particle, p(N), first broadens, indicating a more heterogeneous structure, but then narrows again indicating that the probability to find monomers is reduced to zero, consistent with the observed network structure (Figs 1a,3a). Along path B, that is, increasing c salt (inset), the distributions also shift to larger N and broaden, again indicating the coarsening of the structure. Compared with path A, the p(N) are sharper, indicating that monomers are suppressed at low salt concentrations. Under these conditions, repulsions are long ranged and still substantial at the mean monomer separation. The system avoids this unfavourable situation by forming clusters. Although the repulsion between the clustered monomers is even larger, cluster formation is favoured by the short-ranged attractions and the larger mean cluster separation that reduces cluster-cluster repulsion.
We characterize p(N) by its mean hNi, variance hDN 2 i and maximum N max . With increasing attraction, that is, |V D,min |, they all sharply increase at V min E À 3.0 k B T, where k B T is the thermal energy (Fig. 3b). Furthermore, we consider the distribution p j ð Þ of angles between two successive bonds, j. The distribution is normalized by the solid angle covered by the angle j, that is, such that 2p R p j ð Þsin j dj¼1. It reveals a qualitative change that also occurs at V min E À 3.0 k B T (Fig. 3c). While the peak at j¼p=3, which is related to locally dense tetrahedral packings, grows monotonically with |V min |, a second peak at j¼ 2p=3 develops and becomes pronounced for V min j j43:0 k B T. Furthermore, for these large potential depths, a peak at j¼p emerges, which indicates the formation of straight strings of connected particles. Therefore, based on the mean number of bonds per particle hNi, the different structures on the path towards a gel can be distinguished and quantitatively characterized.
Continuous and DP transitions. We will now show that the particle networks undergo two percolation transitions, a continuous percolation transition and a DP transition. In both cases, percolated networks span the whole system. In continuous percolation, a path along a percolating cluster may contain steps in all directions, including backward steps (Fig. 4a, P). Instead, in the case of DP, only paths along an arbitrarily chosen direction  Table 1). The lower row represents the corresponding states as observed in the simulations. (b) Corresponding pair correlation functions g(r) determined from the three-dimensional data as a function of the particle-particle distance r in units of the particle diameter s, which, for the simulation data, is taken to be the effective particle diameter, that is, s ¼ s eff . Symbols represent experimental data, solid lines simulation data and (only for B1) the light blue solid line simulation data with parameters obtained by an unconstrained fit (for details see Supplementary Note 1). Data have been shifted for clarity. The inset shows the (unshifted) first peaks. Mean number of bonds per particle, hNi, as a function of the minimum of the attractive part of the potential, evaluated at the global minimum, V D,min in units of the thermal energy k B T, and the inverse screening length k normalized by the inverse particle diameter s À 1 , which, for the simulation data, is taken to be the mean particle diameter, that is, s ¼ s 0 , in a semi-logarithmic representation. The background colour represents hNi as obtained by simulations, while the positions and colours of the circles indicate the compositions and hNi, respectively, of the experimental samples, which are labelled. The two white lines indicate two paths, during which either only the polymer concentration ('path A') or only the salt concentration ('path B') are changed. The lower and upper dark solid lines indicate the continuous and directed percolation transitions, respectively, which are taken to occur when the corresponding p(l box ) ¼ 0.2 (Fig. 4c).
are considered whereas possible backward steps in this direction are not considered (Fig. 4a, DP). In the limit of an infinite system size, percolating clusters occur in all directions in the sample. Therefore, on a macroscopic scale no preferred direction exists in a directed percolated network.
The (cumulative) probabilities p(Dx) that a particle participates in paths with a projected length of at least Dx and in directed paths with a length of at least Dx along any arbitrarily chosen direction are shown in the inset and main part of Fig. 4b, respectively, where the solid lines represent simulation data along path A and symbols indicate results from experiments along path B. For both percolation conditions, the probabilities decay to zero in the case of weak attractions (small |V min |). In contrast, for strong attractions plateaus develop, indicating large clusters that span the whole system. The height of the plateaus correspond to the fractions p(l box ) of particles that are part of clusters with the maximum possible length, that is, the box size l box (Fig. 4c). The system size determines l box and also affects the simulations. Nevertheless, within the investigated system sizes only a weak dependence on system size, quantified by the total number of particles M, is observed. In particular, the occurrence of the sharp increases of p(l box ) and their positions are not significantly affected. These sharp increases indicate the continuous and DP transitions.
We determine the continuous and DP transitions in simulations for the complete parameter space of Fig. 2. The locations are chosen to coincide with the percolation probability p(l box ) ¼ 0.2 for systems with 9,856 particles (Fig. 4c). Using this condition, along path A, the continuous percolation transition occurs at V P ¼ À 2.3 k B T and the DP transition is found at V DP ¼ À 3.0 k B T. The latter coincides with the significant changes observed in other parameters; hNi (Fig. 2, dark solid lines) as well as p(N), hDN 2 i, N max and p j ð Þ (Fig. 3). Furthermore, this indicates that the continuous percolation transition occurs at hNiE2 and the DP transition at hNiE3 (Fig. 3b). This suggests a connection between the percolation transitions and hNi. Previous work 32,37 on particles with attractive square-well interactions also observed continuous percolation for hNiE2. This work moreover suggested that rigidity percolation is associated with gelation, which was found to occur at hNiE2.4, that is, just in between our estimates for the onsets of continuous and DP. The possible link between rigidity and DP deserves further investigation.
Furthermore, close to the transitions the values of p(l box ) obey critical power law scaling with V P À V min and V DP À V min , respectively (Fig. 4c,d). Note that up to leading order the probability for a bond between neighbours is proportional to V min and therefore scaling laws for p(l box ) as a function of this probability also hold for p(l box ) as a function of V min . Our data are consistent with theoretical predictions for the critical exponents of the power law scalings, that is, b P ¼ 0.42 for continuous percolation 42,43 and b DP ¼ 0.58 for directed percolation 44 (black lines in Fig. 4c,d; only the intercepts are fitted and the data with almost constant p(l box ) at large V P À V min and V DP À V min , respectively, are not included in the fits). These findings strongly suggest that, on the path towards a gel, continuous and DP networks occur with well-defined transitions to other states as well as between them. Furthermore, this links the formation of gels to a universality class of transitions, namely DP, and hence a well-developed theoretical framework.
Onset of slowdown and ageing in directed percolated systems. The dynamics of the different structures are quantified by the self-intermediate scattering function F and mean squared displacement hDr 2 i (Fig. 5). Diffusive dynamics is observed for individual particles and small clusters, characteristic for fluids. In continuous percolated samples, open network structures are formed (Fig. 1, B2) that are rearranged by local motions requiring only slightly more time of the order of the Brownian time t B (Fig. 5). In contrast, a marked slowdown of the dynamics is observed beyond the directed percolation (DP) transition. DP is characterized by dense clusters with thick strands (Fig. 1, B3) and thus rearrangements are expected to involve long-range, global motions and the breaking of several bonds. Correspondingly, they require long times. The dynamics of the different states hence are consistent with our structural observations. They suggest that the slowdown associated with gelation sets in at the DP transition and the dynamics eventually arrests deeper in the gel region.
The simulation data presented so far have been determined for a fixed waiting time of 300 t B after the initial quench. The effect of the waiting time is illustrated for the dependence of the mean number of bonds per particle, hNi, on the potential depth V min (Fig. 6), which has been linked to the continuous and DP transitions (Fig. 3b). For small attraction strengths, V min j j 3:0 k B T, hNi is independent of the waiting time within the examined time range. However, for V min j j43:0 k B T and thus for DP, hNi is found to increase with waiting time. This indicates that ageing effects are important in directed percolated systems and equilibrium is reached only very slowly. In continuous percolation, in contrast, equilibrium is attained quickly. These findings also indicate the importance of the above-mentioned local and global rearrangement processes in continuous and directed percolated systems, respectively. In addition, it also supports the above conclusion that the DP transition is an equilibrium transition just before the onset of dynamical arrest.
Synaeresis in directed percolated systems. We also investigate the vicinity of walls (Fig. 7). Initially, the particles fill the whole volume, also close to the wall. In fluids and continuously percolated samples this remains like this at least for several weeks (about 10 6 t B ). In contrast, in directed percolated systems, the vicinity of the walls becomes depleted of particles in less than a day (about 3 Â 10 4 t B ), except for very few clusters which remain attached to the cover slip in the experiments, likely due to depletion and van der Waals attractions 45 . The depletion close to the walls indicates that directed connections between particles tend to compact the network. This resembles synaeresis, that is, the macroscopic expulsion of fluid from a gel due to the shrinking of the network. This has been observed in a variety of materials, like gelatin 46 , polysaccharide gels 47,48 , protein gels 49 , organogels 50 , microgels 51 and also weakly attractive colloidal gels 52 . The link between network shrinkage and DP might provide new aspects for the understanding of synaeresis.

Discussion
We investigate gel formation in a system with competing attractive and repulsive interactions. Confocal microscopy experiments on charged colloid-polymer mixtures are combined with Brownian dynamics simulations of particles interacting via the Asakura-Oosawa and Coulomb potentials. Depending on the overall potential minimum V min , which can be varied through the attractive and/or repulsive component of the interactions, different states are identified; fluids of individual particles or clusters as well as continuous and directed percolated networks. The transitions between these states are associated with changes in structural and dynamical parameters, in particular, the number of bonds per particle with significant increases in the mean, variance and most probable value of the distribution of the number of bonds.
The effect of continuous percolation on the dynamics is small. However, DP leads to a significant slowdown of the dynamics. This is attributed to the large number of bonds and concomitant strong confinement of the particles in the attractive potentials of their neighbours, which is also reflected in significant ageing observed in directed percolated systems. It suggests that equilibration is very slow in directed percolated systems, but occurs quickly in continuous percolated systems. Hence, our findings suggest that the slowdown associated with gelation already sets in at the DP transition. Deeper in the gel region, the dynamics then becomes arrested without any additional structural transition occurring.
Our finding that DP precedes the formation of arrested gels, establishes a relation between structural and dynamic features of gels. While previous studies linked gelation in adhesive hard spheres to continuous percolation 28,29,31 , we argue that in colloid-polymer mixtures with competing attractive and repulsive interactions directed rather than continuous percolation indicates the transition to gel states. Moreover, DP is a universality class of transitions. Thus gel formation can be linked to critical behaviour and the corresponding theoretical formalism 53 , which can be exploited to obtain a deeper understanding of gel formation. Moreover, maybe gel collapse under gravity 30 is also initiated by percolation transitions. Our results hence contribute to an improved understanding of the relation between structural and dynamic features of gel-forming systems and hint at their importance for gel collapse and rheological properties.
Furthermore, in directed percolated systems the dense and directed networks appear to lead to a contraction. This results in a detachment from the walls and hence a layer that is depleted of particles appears. This phenomenon resembles synaeresis and might have important implications for applications and industrial products but also, for example, gel collapse and rheological measurements, in particular concerning wall slip 45,[54][55][56] .

Methods
Experiments. The samples consist of spherical PMMA particles, sterically stabilized with poly(12-hydroxystearic acid) polymers and fluorescently labelled with 7nitrobenzo-2-oxa-1, 3-diazole-methyl-methacrylate, as well as linear polystyrene. The particles have a diameter s ¼ 1.72 mm and a polydispersity of 7%, which are determined by static and dynamic light scattering of a very dilute suspension (volume fraction F ¼ 0.005). They are dispersed in a mixture of cis-decalin and cycloheptylbromide, which matches the density and almost matches the refractive index of the particles. We do not observe any indication of anisotropy, sedimentation or creaming of the samples during the observation time (days), indicating a high degree of density matching between solvent and particles. In some samples, creaming is observed after several weeks. In this solvent mixture, PMMA particles acquire a charge 33,57 .
To prepare the colloid stock solution, the particles are sedimented using a centrifuge. The sediment is assumed to be random close packed with a volume fraction F RCP ¼ 0.64 and used to prepare a colloid stock solution with FE0. 40. The polystyrene has a molar mass M w ¼ 3 Â 10 6 kg mol À 1 resulting in a radius of gyration r g ¼ 65 nm (ref. 58) and polydispersity characterized by M w /M n ¼ 1.17. For a dilute mixture, the polymer-colloid size ration 2r g /s ¼ 0.076. Since r g depends on concentration, an effective polymer-colloid size ratio x eff E0.03 is estimated based on the Generalized Free Volume Theory 59-61 (Supplementary Table 1). Different polymer stock solutions are prepared by dispersing dry polymer in the solvent mixture to yield concentrations 4.3 mg ml À 1 rc p ¼ m p r s /m s r18.9 mg ml À 1 with m p the mass of the dry polymer, r s the density and m s the mass of the solvent mixture, respectively. Furthermore, the overlap concentration is calculated using c Ã p ¼3M w =4pN A r 3 g with N A Avogadro's constant. To obtain colloid-polymer mixtures with FE0.2, colloid and polymer stock solutions are mixed with a volume ratio of 50:50 for samples B1-B3, C1, C3 and D1, while dry polymer is added to the diluted colloid stock solution in the case of samples C2, C4 and C5. (The exact composition of each sample is reported in Supplementary Table 1.) To homogenize the samples, the mixtures are left on a vortex mixer for some minutes and then on a flask shaker for at least one day. To vary the screening length k À 1 , different amounts of salt (tetrabutylammoniumchloride, Supplementary Table 1) 33 are added and the mixture gently moved on the flask shaker for three days to help the dissolution of the salt in the organic solvent mixture, which is slow and only very limited. Therefore, unfortunately, the exact amount of dissolved salt is unknown. For a quantitative comparison with the simulations, hence a fitting procedure is applied (Supplementary Note 1).
Within 3 h after mixing, the samples are transferred to a home-built sample cell 62 and imaged with a confocal unit (Visitech VT-Eye and Nikon A1R-MP for the structural and dynamic measurements, respectively) that is mounted on an inverted microscope (Nikon Ti-U and Nikon Ti-E, respectively) equipped with an oil immersion objective (Nikon Plan Apo VC Â 100). In the structural measurements, for each sample, 25-32 image stacks are recorded, each consisting of 151 frames (512 Â 512 pixels) with z steps of 200 nm, corresponding to a volume of 54 Â 54 Â 30 mm 3 , which contains about 7,000 particles. With a rate of 30 slices per second and averaging of 3 slices to obtain 1 frame, a sampling time smaller than 20 s per stack is achieved. For the determination of the dynamics, two-dimensional slices (512 Â 512 pixels) corresponding to an area 63.5 Â 63. are recorded. For each sample, at least two series with 10,000 slices at a rate of 15 slices per second and 4,000 slices at 1 slice per second are taken. The particle coordinates are extracted using standard methods 63 and, where necessary, are refined 62 . The mean squared displacement is calculated from the particle coordinates according to hDr 2 ðtÞi i;t0 ¼h½r i ðt 0 þ tÞ À r i ðt 0 Þ 2 i i;t0 where the average is taken over particles i and starting times t 0 . For the intermediate scattering function Fðk; tÞ¼hexp½ik ÁDr i ðtÞi i;t0 ¼ exp ik Á Dx i þ Dy i þ Dz i ð Þ ½ h i i;t0 the magnitude of the scattering vector is chosen as k ¼ 2pr 1/3 E2.37 mm À 1 . The delay time is scaled with the Brownian time t B ¼s 2 = 3pD 0 ð Þ % 3:1 s for our system with the short-time self-diffusion coefficient D 0 ¼ k B T/g and the friction coefficient of a particle, g. Note that we correct for possible drift in the set-up by subtracting the centre-of-mass motion of all particles in the observation volume.
Computer simulations. We simulate a system consisting of M particles with mean diameter s 0 and a polydispersity of about 7%. The polydispersity is realized by shifting the interaction potentials multiple times corresponding to 17 different particle diameters s i , whose numbers M i are normally distributed with a standard deviation of the diameter of 0.07s and M¼ P i M i . If not stated otherwise, M ¼ 9,856.
Particles interact with a potential V(r ij ) ¼ V C (r ij ) þ V D (r ij ) that is composed of a repulsive screened Coulomb part 64 where s ij ¼ (s i þ s j )/2, and a short-ranged depletion potential according to the continuous Asakura-Oosawa potential 34,65 plus a steep power law potential Motivated by the experiments, x eff ¼ 0.03 is fixed although in our model this does not necessarily imply that the position of the potential minimum is fixed. The hard interaction of particles at contact is approximated by V H ðr ij Þ¼H 0 ðr ij =s ij Þ À 32 (ref. 66). The prefactor W 0 defines the energy scale of the attractive part and the prefactor H 0 the energy scale of the hard-sphere-like potential and is kept constant at H 0 ¼ 0.25 k B T. We cut the potential and the corresponding force at 4s and shift them to zero at the cut-off distance. In our model, the range of the attractive part in the potential may vary while other parameters are changed. The interaction potential exhibits a minimum V min ¼Vðs eff ij Þ with the effective diameter s eff ij slightly smaller than s ij . The effective packing fraction F eff ¼ðp=6VÞ P i M i ðs eff i Þ 3 ¼0:2, with V the volume of the system, is kept constant. For bulk simulations we employ periodic boundary conditions to the cubic simulation boxes. For simulations with parallel flat walls, particles are exposed to an external force, which is defined for the bottom wall as ; zos i =2 0 ; else where we fix F wall s ¼ 10 4 k B T. An analogous force is introduced for the top wall. These forces imply interactions close to hard wall interactions but with a small soft contribution to be suitable for simulations. In our Brownian dynamics simulations particle trajectories are governed by their individual Langevin equations 67 with the friction coefficient g i and hf i t ð Þf T j t 0 ð Þi¼2g i k B Td ij d t À t 0 ð ÞI 3 , where f T j is the transpose of f j , I 3 the three-dimensional unity matrix, g i the friction coefficient of a particle with diameter s i , and d(t À t 0 ) and d ij the Dirac delta and the Kronecker delta, respectively.
As in the experiments, the Brownian time t B ¼ s 2 /(3pD 0 ), with the short-time self-diffusion coefficient D 0 ¼ k B T/g and the friction coefficient of a particle with average diameter, g, is used as time scale. The time steps are 10 À 4 t B or less. Starting from an initially randomly distributed ensemble, the systems are equilibrated for 300 t B before statistics are gathered, unless stated otherwise.
Data availability. The data from experiments as well as the data of the simulations are available from the corresponding author.