A molecular dynamics study of chemical gelation in a patchy particle model
Silvia Corezzi, Cristiano De Michele, Emanuela Zaccarelli, Daniele Fioretto, Francesco Sciortino
aa r X i v : . [ c ond - m a t . s o f t ] F e b A molecular dynamics study of chemical gelation in a patchy particle model
Silvia Corezzi,* ab Cristiano De Michele, cd Emanuela Zaccarelli, cd Daniele Fioretto bc and Francesco Sciortino cd a CNR–INFM Polylab, Universit ` a di Pisa, Largo Pontercorvo 3, I-56127 Pisa, Italy b Dipartimento di Fisica, Universit ` a di Perugia, Via A. Pascoli, I-06100 Perugia, Italy c CNR–INFM CRS Soft, Universit ` a di Roma “La Sapienza”, P. A. Moro 2, I-00185 Roma, Italy d Dipartimento di Fisica, Universit ` a di Roma “La Sapienza”, P. A. Moro 2, I-00185 Roma, Italy (Dated: October 28, 2018)We report event–driven molecular dynamics simulations of the irreversible gelation of hard ellip-soids of revolution containing several associating groups, characterizing how the cluster size distri-bution evolves as a function of the extent of reaction, both below and above the gel point. We findthat in a very large interval of values of the extent of reaction, parameter–free mean–field predictionsare extremely accurate, providing evidence that in this model the Ginzburg zone near the gel point,where non–mean field effects are important, is very limited. We also find that the Flory’s hypothesisfor the post–gelation regime properly describes the connectivity of the clusters even if the long–timelimit of the extent of reaction does not reach the fully reacted state. This study shows that irre-versibly aggregating asymmetric hard–core patchy particles may provide a close realization of themean–field model, for which available theoretical predictions may help control the structure and theconnectivity of the gel state. Besides chemical gels, the model is relevant to network–forming softmaterials like systems with bioselective interactions, functionalized molecules and patchy colloids. I. INTRODUCTION
Irreversible polymerization is a mechanism of self–organization of molecules which proceeds via the forma-tion of covalent bonds between pairs of mutually–reactivegroups.
If monomers with functionality (number f of reactive groups on a monomer) greater than two arepresent, branched molecules grow by reactions and con-vert the system from a fluid of monomers into a wellconnected cross–linked network, giving rise to a chem-ical gelation process. At the gel point, a persistentnetwork spanning the sample first appears; the systemis then prevented from flowing, yet not arrested on amesoscopic length scale. The development of a networkstructure results, for example, from step polymerization,chain addition polymerization and cross–linking of poly-mer chains. The same phenomenon is also observed incolloids and other soft materials when the thermodynam-ics and the molecular architecture favor the formation ofa limited number of strong interactions (i.e., with at-traction strength much larger than the thermal energy)between different particles. Chemical gelation has beenextensively studied in the past, starting from the pio-neering work of Flory and Stockmayer who developedthe first mean–field description of gelation, providing ex-pressions for the cluster size distribution as a functionof the extent of reaction and the critical behavior of theconnectivity properties close to gelation. More appro-priate descriptions based on geometric percolation con-cepts have, in the late seventies, focused on the non–meanfield character of the transition, which reveals itself nearthe gel point, extending to percolation the ideas devel-oped in the study of the properties of systems close to asecond–order critical point. Several important numericalstudies, —most of them basedon simulations on lattice —have focused on the criticalbehavior close to the percolation point, providing evi-
FIG. 1: Graphic description of the A and B particles (left)and snapshot of the simulated system (right). The centers ofthe small spheres locate the bonding sites on the surface ofthe hard–core particle. dence of the percolative nature of the transition and accu-rate estimates of the percolation critical exponents. As incritical phenomena, a crossover from mean–field to per-colation behavior is expected close to the gel transition. But, how the microscopic properties of the system con-trol the location of the crossover (i.e., how wide is the re-gion where the mean–field description applies) and howaccurate is the mean–field description far from the per-colation point is not completely understood. Anotherimportant open question regards the connectivity prop-erties of chemical gels well beyond percolation. Evenin the mean–field approximation, several possibilities forthe post–gel solutions have been proposed, based on dif-ferent assumptions on the reactivity of sites located onthe infinite cluster.
Different propositions predict dif-ferent cluster–size distributions above the gel point anda different evolution with time for the extent of reaction.Here we introduce a model inspired by step-wise polymerization of bifunctional diglycidyl–ether ofbisphenol–A ( B particles in the following) with penta-functional diethylenetriamine ( A particles). To incor-porate excluded volume and shape effects, each type ofmolecule is represented as hard homogeneous ellipsoidof appropriate length, whose surface is decorated in apredefined geometry by f identical reactive sites per par-ticle (see Figure 1). In this respect, the model is alsorepresentative of colloidal particles functionalized with alimited number of patchy attractive sites, where theselectivity of the interaction is often achieved buildingon biological specificity. The off–lattice evolutionof the system is studied via event–driven molecular dy-namics simulations, using a novel code which specificallyextends to ellipsoidal particles the algorithm previouslydesigned for patchy spheres. Differently from previousstudies, we do not focus on the critical properties closeto the gel–point but study in detail the development ofthe irreversible gelation process and the properties of thecluster size distribution in the pre– and post–gelationregime.We find that the dynamic evolution of the system pro-duces an irreversible (chemical) gelation process whoseconnectivity properties can be described, in a verylarge window of the extent of reaction, with the Flory–Stockmayer (FS) predictions.
This offers to us thepossibility to address, in a well controlled model, the ki-netics of the aggregation and to evaluate the extent ofreaction at which the breakdown of the Flory post–gelsolution takes place.
II. METHOD
We study a 5:2 binary mixture composed of N A = 480ellipsoids of type A and N B = 1200 ellipsoids of type B ,for a total of N = 1680 particles. A particles are modeledas hard ellipsoids of revolution with axes a = b = 2 σ and c = 10 σ and mass m ; B particles have axes a = b = 4 σ and c = 20 σ , mass 3 . m . Simulations are performed ata fixed packing fraction φ = 0 .
3. Five (two) sites arerigidly anchored on the surface of the A ( B ) particles,as described in Fig. 1. Sites on A particles can only re-act with sites on B particles. Every time, during thedynamic evolution, the distance between two mutually–reactive sites becomes smaller than a predefined distance δ = 0 . σ , a new bond is formed between the particles.To model irreversible gelation, once a bond is formed,it is made irreversible by switching on an infinite bar-rier at distance r ijAB = δ between the sites i and j in-volved, which prevents both the formation of new bondsin the same sites and the breaking of the existing one.Hence, the newly formed bond cannot break any longer,and the maximum distance between the two reacted sitesis constrained to remain smaller than δ . Similarly, the FIG. 2: Time dependence of the fraction of bonds p . Sym-bols: simulation results (averaged over 40 independent real-izations). For the chosen stoichiometry, p coincides with thereacted fraction of A reactive sites, i.e. the A conversion, orequivalently with the reacted fraction of B sites, i.e. the B conversion. p = 1 would indicate that all possible bondingsites have reacted. Time is measured in arbitrary units. Line: p ( t ) = kt/ (1 + kt ), with the fit–parameter k fixing the timescale. This functional form is expected when any pair of re-active groups in the system is allowed to react, but loops donot occur in finite size clusters. two reacted sites cannot form further bonds with avail-able unreacted sites. The composition of the system andthe particle functionality are such that the reactive sitesof type A and B are initially present in equal number, f A N A = f B N B , which in principle allows the formationof a fully bonded state in which all the sites have reacted.This offers a way to properly define the extent of reac-tion as the ratio p between the number of bonds presentin a configuration and the maximum number of possiblebonds f A N A .Between bond–formation events, the system propa-gates according to Newtonian dynamics at temperature T = 1 .
0. As in standard event–driven codes, the config-uration of the system is propagated from one collisionalevent to the next one. Note that temperature only con-trols the time scale of exploration of space, by modulatingthe average particle’s velocity. An average over 40 inde-pendent starting configurations is performed to improvestatistics.
III. RESULTS
In the starting configurations no bonds are present byconstruction. As a function of time, the fraction p offormed bonds —a measure of the state of advancement ofthe reaction— increases monotonically, until most of theparticles are connected in one single cluster (Figure 2).As a result, p saturates around 0 .
86, despite the factthat an equal number of reactive sites of type A and B is initially present in the system.Flory and Stockmayer laid out the basic relations be-tween extent of reaction and resulting structure in steppolymerizations, on the assumptions that all functionalgroups of a given type are equally reactive, all groupsreact independently of one another, and that ring for-mation does not occur in molecular species of finite size.Only when p exceeds a critical value p c infinitely largemolecules can grow. In this respect the FS theory de-scribes the gelation transition as the random percolationof permanent bonds on a loopless lattice. The presentmodel satisfies the conditions of equal and independentreactivity of all reactive sites. The absence of closedbonding loops in finite size clusters is not a priori im-plemented; as we will show in the following, however,such a condition —favored by the poor flexibility of thebonded particles and their elongated shape, the absenceof an underlying lattice and the asymmetric location ofthe reactive sites— is valid in a surprisingly wide regionof p values.The FS theory predicts the p dependence of the clus-ter size distribution in the very general case of a mixtureof monomers bearing mutually reactive groups. In thepresent case, the number n lm of clusters containing l bi-functional particles and m pentafunctional ones can bewritten as n lm = N B N A p l + m − (1 − p ) m +2 w lm (1) w lm = (4 m )!( l − m + 1)!(4 m − l + 1)! m !and the number of clusters of size s is obtained by sum-ming over all contributions such that l + m = s , i.e., n s = P lm,l + m = s n lm . As shown in Figure 3a on increas-ing p , the n s distribution becomes broader and broaderand develops a power–law tail. The theory predicts agelation transition when p c = 1 / p ( f A − f B −
1) =0 . Even close to p = 0 .
5, the FS prediction —whichconforms to the prediction of random percolation on aBethe (loopless) lattice where n s ∼ s − . at the perco-lation threshold— is consistent with the numerical data.On further increasing p (Figure 3b), the distribution offinite size clusters progressively shrinks, and only smallclusters survive. Data show that Eq 1, with no fitting pa-rameters, predicts rather well the numerical distributionat any extent of polymerization, both below and abovethe point where the system is expected to percolate, in-cluding details such that the local minimum at s = 2.To compare with the mean–field prediction of gelationat p c = 0 .
5, we examine the connectivity properties ofthe aggregates for each studied value of p , searching forthe presence of clusters which are infinite under peri-odic boundary conditions. We find that configurationsat p = 0 . ± .
008 have not yet developed a percolat-ing structure while configurations at p = 0 . ± . p c = 0 . ± . FIG. 3: Distribution of finite size clusters n s for different frac-tion of bonds p (a) below and (b) above percolation. Pointsare simulation data and lines are the corresponding theoret-ical curves according to FS. The dashed line represents thepower law n s ∼ s − . . to the infinite (percolating) network N ∞ constitutes the gel , while the soluble material formed by the finite clus-ters which remain interspersed within the giant networkconstitutes the sol . Figure 4a shows that the fraction ofgel P ∞ = N ∞ /N and even its partition between particlesof type A ( P A, ∞ = N A, ∞ /N ) and B ( P B, ∞ = N B, ∞ /N )calculated according to the FS theory, properly rep-resent the simulation results throughout the polymer-ization process. Indeed, the proportion of B particlesto A particles in gel and in sol is a function of p (seeinset). The relative amount of B particles in the sol( N B,sol /N A,sol ) increases as a consequence of the prefer-ential transfer of the A particles (having more reactivesites) to the gel, in a way that the fraction p sol of sites B in the sol that have reacted (extent of reaction in the sol)differs from the total fraction p of sites B reacted (extentof reaction in the system). The constitution of the sol(Figure 3(b)) results to be the same as that of a smallersystem made of N A,sol particles of type A and N B,sol particles of type B reacted up to the extent p sol . The evolution of the cluster size distribution can bequantified by the number ( x n ) and weight average ( x w )cluster sizes of the sol, defined as x n = P s sn s / P s n s and x w = P s s n s / P s sn s . The numerical results and FIG. 4: (a) Gel fraction P ∞ and its partition between particlesof type A ( P A, ∞ ) and B ( P B, ∞ ) vs the fraction of bonds p (i.e. the extent of reaction in the system). The inset shows theproportion of B particles to A particles in gel ( N B, ∞ /N A, ∞ — left axis) and in sol ( N B,sol /N A,sol — right axis) vs p . (b)Number and weight average cluster size ( x n and x w ) prior togelation and for the sol after gelation vs the fraction of bonds p . (c) Relation between the number of finite size clusters(molecules in the sol) n sol and the fraction of bonds p . Theinset shows the number of loops n loop vs p . In all panels,symbols are simulation results and solid lines FS predictions. the FS theoretical predictions are shown in Figure 4b.Both averages increase before gelation; then, they regressin the sol existing beyond the gel point, since large clus-ters are preferentially incorporated into the gel network.While x n increases only slightly up to the gel point, neverexceeding 3.5, x w increases sharply in proximity of p c as well as sharply decreases beyond this point, consis-tently with the fact that x w is singular at percolationbeing dominated by large clusters. Again, simulationdata agree very well with FS predictions. Discrepan-cies between theory and simulation —which reveal themean–field character of the FS theory— only concern therange of p very near p c , suggesting that for this modelthe crossover from mean–field to percolation is very closeto the gel point — i.e., the Ginzburg zone near thegel point, where non–mean field effects are important, isvery limited. A finite–size study very close to the criticalpoint would be requested to accurately locate the per- colation point and the critical exponents, a calculationbeyond the scope of the present work.From a physical point of view, the change from mean–field to percolation universality class is rooted in the pres-ence of bonding loops in the clusters of finite size, whichpre-empts the possibility to predict the cluster size distri-bution. The realistic estimate of the percolation thresh-old and the agreement between theory and simulation(Fig. 3) suggest that the present model strongly disfa-vors the formation of loops in finite clusters, at leastfor cluster sizes probed in this finite–size system. As atest, we evaluate the total number of finite (sol) clusters n sol = P s n s as a function of the extent of reaction. Iffinite clusters do not contain closed loops, n sol equals thenumber of particles in the sol minus the number of bonds,since each added bond decreases the number of clustersby one. This applies equally to the system preceding gela-tion, or to the sol existing beyond the gel point. Thus, at p < p c (pre–gelation) the relation between n sol and p islinear, i.e. n sol = N − N B p . At p > p c (post–gelation), n sol can be calculated as n sol = N sol − N B,sol p sol , where N sol is the number of particles in the sol fraction ( N B,sol of which bear reactive sites of type B ), and p sol = p is thereacted fraction of sites B in the sol. Hence, the relationbetween n sol and p crosses to a nonlinear behavior, sothat the number of clusters becomes one when p = 1. Asshown in Figure 4c, the number of finite clusters found inthe simulation data conforms to the theoretical expecta-tion for all p values, both below and above the gel point.Hence, as a first approximation, loops are only present inthe infinite (percolating) cluster and do not significantlyalter the distribution of the finite size clusters, both be-low and above percolation. The difference between n sol found in simulation and the value predicted by the FStheory counts the number of loops in the sol, n loop . Sucha quantity is shown in the inset of Figure 4c. The maxi-mum value of n loop , achieved for p ∼ p c , corresponds to0 .
2% of the total number of bonds. This demonstratesthat intramolecular bonds within finite clusters can beneglected, consistent with the Flory hypothesis for thepost–gelation regime . Figure 4c also shows that thelinear relation between n sol and p is valid also after thegel point (up to p ≈ . on the polymer-ization of bifunctional diglycidyl–ether of bisphenol–Awith pentafunctional diethylenetriamine, also suggestingthat the number of cyclic connections in the infinite clus-ter is negligible well above p c .As a further confirmation of the absence of closed loopswe compare the time evolution of p with the predictionof the mean–field kinetic modeling of polymerization,based on the solution of the Smoluchowski coagulationequation. For loopless aggregation, p ( t ) is predictedto follow p ( t ) = kt kt , (2)where the fit–parameter k , which has the meaning of abond kinetic constant, fixes the time scale of the aggre-gation process. The time evolution of p is found to per-fectly agree with the theoretical predictions (see Fig-ure 2) up to p ≈ .
6, i.e. beyond p c . While the predictionwould suggest that p ( t → ∞ ) = 1 (dash line in Figure 2),the simulation shows that the formation of a percolatingstructure prevents the possibility of completing the chem-ical reaction, leaving a finite number of unreacted sitesfrozen in the structure. As shown above (Figure 3), evenin this frozen state the cluster size distribution is pro-vided by the Flory’s post–gel hypothesis. Such a featureis not captured by the mean–field Smoluchowski equationin which spatial information in the kernels are neglected. IV. CONCLUSIONS
A binary mixture of patchy hard ellipsoids undergoingchemical gelation displays a very large interval of the ex-tent of reaction in which parameter–free mean–field pre-dictions are extremely accurate. The connectivity prop-erties of the model are properly described — withoutany fitting parameter — both below and above perco-lation by the mean–field loopless classical FS theory.
The mean–field cluster size distribution for the sol com-ponent is found to be valid for all values of the extent ofreaction, both below and above the gel point, suggestingthat for the present model, the Flory’s hypothesis for thepost–gelation regime properly describes the irreversibleaggregation phenomenon, despite the explicit considera-tion of the excluded volume.The absence of loops in finite size clusters, which isnot assumed by the model, results from the specific ge-ometry of the bonding pattern and by the presence ofthe excluded volume interactions, disfavoring the forma-tion of ordered bonding domains. Hence, the geometryof the particles and the location of the reactive sites onthem may play a significant role in the stabilization ofthe mean–field universality class with respect to the per-colation universality class, locating the crossover be-tween the two classes very close to the gel point. Thepresent study shows that irreversibly aggregating asym-metric hard–core patchy particles, even if excluded vol-ume effects are properly taken into account, may providea close realization of the FS predictions in a wide rangeof p values. The model thus offers a starting point —forwhich theoretical predictions are available— for further investigations of the gelation process and for a more pre-cise control over the structure and connectivity of thegel state. In particular, a full and detailed structural in-formation can be known along with the dynamics of thesystem, which is potentially useful to investigate the rela-tion between structural heterogeneity and heterogeneousdynamics, and to shed light on the microscopic aspectsof the dynamic crossover from short to long relaxationtimes, during irreversible polymerization.While the structural properties are all well–describedby the FS theory, the evolution of the extent of reaction,modeled via the coagulation Smoluchowski equation, isproperly described by the theory only in the pre–gelationregion. After gelation, kinetic constraints due to the ab-sence of mobility of the reactive sites anchored to thepercolating cluster or to smaller clusters trapped insidethe percolating matrix prevent the completion of the re-action and the extent of reaction freezes (to p ≈ . p value.In the present model, the entire polymerization processproceeds via a sequence of FS cluster size distributions,determined by p ( t ). Recently, it has been shown thatthe FS theory properly describes also equilibrium clus-tering in patchy particle systems when p is a function oftemperature and density. It is thus tempting to spec-ulate that for loopless models, irreversible evolution canbe put in correspondence with a sequence of equilibrium states which could be sampled in the same system for fi-nite values of the ratio between temperature and bondingdepth. If this is indeed the case, chemical gelation couldbe formally described as a deep quench limit of phys-ical gelation. This correspondence would facilitate thetransfer of knowledge from recent studies of equilibriumgels to chemical ones. Concepts developed for irre-versible aggregation of colloidal particles, like diffusion–and reaction–limited cluster–cluster aggregation, couldbe connected to chemical gelation. Work in this direc-tion is ongoing.We acknowledge support from MIUR-PRIN. We thankP. Tartaglia for interesting discussions. P. J. Flory,
Principles of Polymer Chemistry , Cornell Uni-versity Press, London, 1953. M. Rubinstein and R. H. Colby,
Polymer Physics , OxfordUniversity Press Inc., New York, 2003. W. Burchard,
Polym. Bull. , 2007, , 3-14. R. J. Young and P. A. Lovell,
Introduction to polymers ,Chapman and Hall, New York, 1991. J. P. Pascault, H. Sauterau, J. Verdu and R. J. J. Williams,
Thermosetting Polymers , Marcel Dekker, New York, 2002. W. H. Stockmayer,
J. Polym. Sci. , 1952, , 69. D. Stauffer,
J. Chem. Soc., Farad. Trans. II , 1976, ,1354. P. Manneville and L. de Seze, in
Numerical Methods inthe Study of Critical Phenomena , ed. J. Della Dora, J.Demongeot and B. Lacolle, Springer-Verlag, Berlin, 1981. H. J. Herrmann, D. P. Landau and D. Stauffer,
Phys. Rev.
Lett. , 1982, , 412. R. B. Pandey and D. Stauffer,
Phys. Lett. , 1983, , 511;Y. Liu and R. B. Pandey,
J. Phys. II France , 1994, , 865-872; Y. Liu and R. B. Pandey, Phys. Rev. E , 1996, ,6609; Y. Liu and R. B. Pandey, Phys. Rev. B , 1997, ,8257-8266. J. P. Clerc, G. Giraud, J. Roussenq, R. Blanc, J. Carton,E. Guyon, H. Ottavi and D. Stauffer,
Ann. Phys. , 1983, ,1. R. Bansil, H. J. Herrmann and D. Stauffer,
Macro-molecules , 1984, , 998. Y. K. Leung and B. E. Eichinger,
J. Chem. Phys. , 1984, , 3887. A. M. Gupta, R. C.Hendrickson and C. W. Macosko,
J.Chem. Phys. , 1991, , 2097. D. Lairez, D. Durand and J. R. Emery,
J. Phys. II France ,1991, , 977-993. J. C. Gimel, T. Nicolai, D. Durand and J. M. Teuler,
Eur.Phys. J. B , 1999, , 91; J. C. Gimel, T. Nicolai and D.Durand, J. Sol-Gel Sci. Tech. , 1999, , 129. D. Vernon, M. Plischke and B. Jo´os,
Phys. Rev. E , 2001, , 031505. E. Del Gado, L. de Arcangelis and A. Coniglio,
Phys. Rev.E , 2002, , 041803. C. P. Lusignan, T. H. Mourey, J. C. Wilson and R. H.Colby,
Phys. Rev. E , 1999, , 5657. M. Rubinstein and A. V. Dobrynin,
Curr. Opin. Coll. Int.Sci. , 1999, , 83. P. G. J. van Dongen,
J. Stat. Phys. , 1997, , 1273. S. Corezzi, D. Fioretto and J. M. Kenny,
Phys. Rev. Lett. ,2005, , 065702. S. C. Glotzer and M. J. Solomon,
Nature Materials , 2007, , 557-562. A. L. Hiddessen, S. D. Rotgers, D. A. Weitz and D. A. Hammer,
Langmuir , 2000, , 9744; A. L. Hiddessen, D.A. Weitz and D. A. Hammer, Langmuir , 2004, , 71. C. A. Mirkin, R. L. Letsinger, R. C. Mucic and J. J.Storhoff,
Nature , 1996, , 607-609. T. Schmatko, B. Bozorgui, N. Geerts, D. Frenkel, E. Eiserand W. C. K. Poon,
Soft Matter , 2007, , 1-5. C. De Michele, S. Gabrielli, P. Tartaglia and F. Sciortino,
J. Phys. Chem. B , 2006, , 8064-8079. D. Stauffer and A. Aharony,
Introduction to percolationtheory , Taylor & Francis, London, 1992. D. R. Miller and C. W. Macosko,
Macromolecules , 1976, , 206. D. R. Miller, E. M. Valles and C. W. Macosko,
Polym.Eng. Sci. , 1979, , 272. S. Corezzi, L. Palmieri, J. M. Kenny and D. Fioretto,
J.Phys.: Condens. Matter , 2005, , S3557. R. Volponi, S. Corezzi and D. Fioretto,
Macromolecules ,2007, , 3450. R. M. Ziff, J. Stat. Phys. , 241 (1980). H. Galina, and J. B. Lechowicz, Adv. Polym. Sci. , 135(1998) D. Stauffer, A. Coniglio and M. Adam,
Adv. Polym. Sci. ,1982, , 193. S. Corezzi, L. Comez, D. Fioretto, G. Monaco and R. Ver-beni,
Phys. Rev. Lett. , 2006, , 255702. S. Corezzi, D. Fioretto and P. A. Rolla,
Nature , 2002, , 653-656; S. Corezzi, D. Fioretto, D. Puglia and J.M. Kenny,
Macromolecules , 2003, , 5271-5278. E. Bianchi, P. Tartaglia, E. La Nave and F. Sciortino,
J.Phys. Chem. B , 2007, , 11765-11769. E. Bianchi, J. Largo, P. Tartaglia, E. Zaccarelli and F.Sciortino,
Phys. Rev. Lett. , 2006, , 168301. E. Zaccarelli,
J. Phys.: Condens. Matter , 2007,19