Distribution Cutoff for Clusters near the Gel Point.

The mechanical and dynamic properties of developing networks near the gel point are susceptible to the distribution of clusters coexisting with percolating networks. The distribution of cluster numbers follows a broad power law, wrapped by a cutoff function that rapidly decays at a characteristic size. The form of the cutoff function has been speculated based on known results from lattice percolation and, in certain cases, solved. We obtained this cutoff function from simulated dynamic clusters of polymeric precursor chains using a hybrid Monte Carlo algorithm. The results obtained from three different precursor chain lengths are consistent with each other and are consistent with the expectation from lattice percolation.

Introduction

Much of the recent interest in dynamic networks was spurred by the need to develop healable and reprocessable elastomers and by the study of associative systems containing long-lived ionic bonds such as ionomers.[1−3] To design materials with high mechanical strength, the cross-linking density needs to be high, whereas to enhance stretchability, the cross-linking density needs to be low (Figure ). Materials with low modulus are sometimes advantageous, for instance, in soft electronics applications.[4] As a result, a balance of strength and stretchability is needed, which often places materials in the gelation window, where the variation of these two properties is the greatest. Elucidating the molecular design rules that correlate backbone stiffness, monomer bulkiness, molecular weight of the precursor chain, and the lifetime of physical cross-links with the mechanical and physical properties of the dynamic network is thus instructional to the design and engineering of soft elastomers that are healable or reprocessable.[5,6]

Figure 1

Variation of modulus,[16] cluster lifetime,[12,13] and stretchability (following ref (2) while assuming that the backbones of the clusters are force-bearing) in the gelation window and above the gel point, and schematics for tree-like clusters (s = 1, 2, 3) and a defective cluster (s = 7) containing a self-loop. The stretchability is estimated as the end-to-end distance of a strand between two cross-links normalized by the unperturbed size.[2] The mean-field to critical regime crossover is ϵG ∝ N–1/3,[20] and the effective breakup sets in for |ϵ| < ϵc.[12,13]

Molecular theories for dynamic networks[7−17] have been developed within the same framework as for permanent networks,[18] which is rooted in percolation theory.[19] The sticky Rouse[8,16,17] and sticky reptation[9] models are commonly used for studying the dynamic behaviors of these systems, and the near-critical gels have been extensively studied.[10,11,14,15] The main results concerning the regime close to the gel point are that the structural and dynamic properties can be captured by a set of closely related scaling exponents for the degree of gelation ϵ. Let the probability of a monomer forming a cross-link be p, and the (mean-field) value of p at the gel point be , where N is the number of monomers per precursor chain.[18] The degree of gelation is defined as . The condition ϵ = −1 corresponds to a melt of precursor chains, while ϵ = 1 corresponds to the completion of gelation.

Within the gelation window defined by −1 ≤ ϵ ≤ 1, which is of primary interest for this work, the system contains a mix of precursor chains, small clusters formed by linked precursor chains and, when ϵ > 0, a fraction of chains belonging to the percolating network. The density of clusters containing s precursor chains obeys a two-parameter scaling ansatz:[18,19]The cutoff function f is constant when s/s* is small and drops to zero otherwise. The exponent τ captures the power-law distribution for cluster numbers, and σ characterizes the typical cluster size formed at a given ϵ. The cluster distribution is the basis for understanding the structural and dynamic features[18] in the vulcanization process studied by Flory[21,22] and Stockmayer.[23] These concepts and analyses can be applied to other associative systems, including colloidal particles,[24] supramolecular clusters,[25,26] and polymer-grafted nanoparticles.[27] Furthermore, as a particular example of a percolation class, the detailed cluster distribution functions obtained from polymeric systems allow generic features from percolation theory to be tested.[19,28]

For gelation in three dimensions (d = 3), τ = 5/2 and σ = 1/2 in the mean-field gelation regime (|ϵ| > N–1/3), and τ = 2.18 and σ = 0.45 in the critical gelation regime (|ϵ| < N–1/3).[18,20] The mean-field regime is obtained from the standard Flory–Stockmayer theory,[18] and the critical regime emerges because the cluster size grows rapidly as the gel point is approached, which eventually overfills space. Consequently, beyond a threshold length, the clusters become space-filling, following the critical exponents.[18] The crossover N–1/3 between these two regimes was first predicted by de Gennes.[20]

The exponents τ and σ determine all other static exponents.[18,19] For instance, the gel fraction is given by Pgel = ϵβ (ϵ > 0), with . The size of clusters is described by the fractal dimension D which, using the space-filling argument,[19] can be shown to be given by in the mean field regime, and by in the critical regime. The correlation length ξ of the system is the size of the largest cluster s*, which scales as , with . The pervaded volume of a cluster of s* precursor chains is ξ3.

The transition from mean-field to critical scaling happens at the Ginzburg number ϵG.[18] Setting the number of clusters of s* precursor chains that can cohabit its pervaded volume to one gives , where is the invariant degree of polymerization. The dependence on N is identical to de Gennes’s original result[20] (Figure ). The dependence on the monomer volume Ω and statistical segmental length b determines how monomer bulkiness and chain stiffness affect the crossover between the mean field (|ϵ| > ϵG) and the critical (|ϵ| < ϵG) regimes. One way to see the effect of stiffness is by replacing the Gaussian statistics with worm-like chain statistics. Setting the number of overlapping strands to one gives , where NK is the number of Kuhn strands per precursor chain. This expression reduces to for flexible chains, whereas it becomes for stiff chains.

Although much is known of the scaling exponents τ and σ for the cluster distributions n(s), the understanding of the cutoff function f(s/s*) remains limited.[18] It is known to exhibit a peak below the gel point, and to decay rapidly further away from the gel point on both sides. Except in certain limiting cases,[18] only numerical results are reported for lattice percolation. For gelation of linear precursor chains, the shape of the cutoff function appears to be absent in the literature.

The aim of the current work is two-fold. First, we present and benchmark a coarse-grained model for the gelation of linear precursor chains. Second, based on the cluster number distributions sampled for multiple precursor chain lengths and multiple boxes, we obtain the cutoff function and examine whether or not a single form applies to different chain lengths. The next section explains the details of our simulation model, including the cross-linking protocol. The following section reports the main results, including the verification of equilibration of the precursor chains, the calibration of the cross-linking rate, the identification of the gel point and, finally, the discussion of the cutoff function. The last section summarizes our main findings.

Model

We use the bead–spring model for precursors.[29] The nonbonded beads interact via the truncated, purely repulsive Lennard-Jones (LJ) potential, , for 0 < r < 21/6σ. Here, the strength of the interaction is set to ε = kBT, and σ is the range of the LJ interaction, which is the default length unit. The bonded beads first interact via the FENE bond potential,[29], where R0 = 1.5 is the maximum bond length and K = 30 is the bonding strength, in addition to an LJ potential. We used cubic simulation boxes with periodic boundary conditions imposed along all three directions. The number of beads per chain is N = 12, 25, or 50, and the number of chains per simulation box is denoted M. Various box volumes are simulated to investigate the finite size effects. The box volume V is set to multiples of V0, with V0 = 14.313. The bead number density ρ = NM/V is set to 0.85 for all systems.

Before the cross-linking steps, the precursor chains are equilibrated following the approach outlined in the literature.[30] The initial configurations are generated by randomly placing the beads. A soft dissipative particle dynamics (DPD) potential is applied for all bead pairs to allow for rapid relaxation. The strength of the DPD potential is then gradually increased until it is comparable to the LJ potential. Next, the DPD potential is replaced with the LJ potential, and a long molecular dynamics (MD) simulation is performed for an isothermal re-equilibration under a Nosé-Hoover thermostat. Finally, the distributions of the end-to-end distance and the radius of gyration of precursor chains are monitored to confirm that the expected random walk statistics is recovered.

Configuration-Biased Cross-Linking

Reversible cross-links or bonds are introduced to the equilibrated melt as follows. Every bead in the system can form a reversible bond with all other beads in its neighborhood. The bond potential is the same as the FENE potential for intramolecular bonds. The number of such bonds is mediated by the chemical potential of the reversible bond μ. At each MC step, the simulation chooses to create or remove a bond with probabilities Pf and Pr = 1 – Pf, respectively. The probability of bond creation or removal is designed according to the Metropolis rule. To facilitate the discussion, we denote the state with one bond removed as the configuration {o}, and the state with one bond formed as the configuration {n}. The MC algorithm moves from {o} → {n} or {n} → {o} are designed by adapting the configuration-biased algorithm,[31] in order to enhance the success rate. The probability of acceptance is designed to fulfill the detailed balance condition.

The bond-creation move begins with randomly selecting a monomer among the monomers capable of forming cross-links, which is denoted i (Figure ). Then all nc neighboring monomers within a cutoff distance Rc = 1.5 are identified, one of which is selected for a bonding attempt. The selection is based on the Boltzmann weight; thus the probability that the bead j is selected is , where r is the distance between beads i and j, u(r) is the bonding energy given by the FENE potential, and is the sum of Boltzmann weights for all the neighbors, i.e., the Rosenbluth factor.[31] The combination of these steps gives the probability of selecting a potential bond that links bead i to bead j as .

Figure 2

Probability flow in the bond-generation algorithm. The number of beads that are capable of forming new bonds is . The distance cutoff for selecting bonding partners is Rc = 1.5. The weighting factors W and W are the sums of the Boltzmann factors for the bonds that can be potentially generated by starting from either bead i or bead j.

The same bond could have been created by first selecting bead j, then identifying bonding partner bead i from its close neighbors (Figure ), which leads to the probability , with W defined from the close neighbors of bead j, in full analogy with W.

To fix the acceptance rate of bond creation while satisfying the detailed balance condition, we have to consider also bond removal which, as mentioned above, is attempted with probability Pr. Let the number of dynamic bonds in configuration {o} be nb. Then the configuration {n} has nb + 1 dynamic bonds, among which the bond between beads i and j is selected for removal with probability . The probability of selecting a given bond to remove during the reverse move is .

Collecting the probabilities for the two bond creation paths (Figure ) and the bond removal step, we have the following detailed balance condition:in which Peq(o) and Peq(n) are the equilibrium probabilities for configurations {o} and {n}, given by the respective Boltzmann weights. Their ratio is solely governed by the bonding energy, i.e., , where z = e is the activity of the dynamic bond. Then the probability of the bond creation move is given byand that for the bond elimination move can be written , in which is an average of the Rosenbluth factors along the two bond generation paths. For this study, we set the attempt frequency ratio to Pf = Pr = 0.5.

The Monte Carlo moves described above are incorporated into an MC/MD algorithm, with each iteration of the algorithm consisting of 1000 MD timesteps followed by 200 attempted MC moves. The MC/MD algorithm is implemented in LAMMPS.[32]

Results

Equilibration of Precursor Chains

To ensure that the melts of precursor chains have been equilibrated, we examined the distributions of the end-to-end vector R and the radius of gyration . The results are shown in Figure and are compared to the predictions of random walk statistics. The prediction for the end-to-end distance R ≡ |R| is given by , in which b is the statistical segmental length. The prediction for the radius of gyration has been derived by Fixman.[33] In both cases, when the argument R or is normalized by N1/2b or Nb2/6, the distributions are expected to be universal.

Figure 3

Distributions of end-to-end distance P(R) and radius of gyration for N = 6, 12, 25, 50, 100. The red curves are predictions using the random walk statistics. The Fixman expression was used for . With increasing N, the distributions converge to the theoretical predictions. Insets: scaling of averages and with molecular weights.

The simulated distributions P(R) and are not identical to theoretical predictions, but the difference decreases as N increases (N = 6, 12, 25), until eventually converging toward the theoretical curves for N = 50, 100. This trend is expected, as the shorter chains are more susceptible to liquid-state packing, which is vividly shown in the spiky patterns of the P(R) data for N = 6. However, these discreteness effects are mild as, even for N = 12 and 25, the data trace the theoretical expectation closely. The random walk scaling is further confirmed in the inset of Figure , whereby the anticipated linear-dependence of the averages and on N is confirmed.

The statistical segmental length b is normally extracted from the scaling of R2 or with N. Because of the modest molecular weights used, we considered the correction to the Gaussian statistics in our estimates. Recent work[34−37] has revisited the conformational statistics of polymers in dense melts, and showed that the deviation from the random walk statistics is dominated by the effects of the correlation hole.[38] Because a polymer spans a finite range, its immediate neighborhood of order always contains certain amount of monomers on the same chain, whose concentration scales as . As a result, the differences between the simulated distributions P(R) and and the random walk predictions are of order N–1/2. This is indeed confirmed by the collapse of the difference normalized by (N – 1)1/2, shown in Figure . These discrepancies modify the scaling of the averages and as well. It has been shown that the N dependence fits to the form .[34−36] In full analogy, is described by . The accuracy of these improved forms is confirmed in the inset of Figure . From the intercept of a straight line fitting, we inferred that the statistical segmental length is b = 1.32. This value was used to normalize the simulated R and in Figures and 4.

Figure 4

Normalized difference between the simulated distributions P(R) and , and the predictions from the random walk statistics for N = 12 (blue), 25 (green), 50 (yellow), and 100 (orange). Insets: (left) and (right) against for N = 6, 12, 25, 50, 100.

Cross-Linking Rate of Dynamic Bonds

The cross-linking probability of dynamic bonds controls the degree of gelation. In our model, the number of dynamic bonds is controlled by the chemical potential μ. For each μ value, the probability p that a monomer forms a dynamic bond is calculated fromwhere the total monomer number is the product NM. Although our model allows monomers to form multiple dynamic bonds, the average number of bonds formed by a cross-linked monomer is close to one, so eq accurately describes the bonding rate (the distribution of cross-links along the chain and the average number of cross-links per bead are shown in Figure S1 in the Supporting Information). The variation of the cross-linking probability with μ gives the binding curve p(μ) shown in Figure . Over a wide range of chemical potential, the binding curve is found to be independent of the molecular weight and the system size, indicating that the intrinsic local cross-linking equilibrium is captured by the MC moves. The dashed line shows a simple quadratic fitting that allows us to target the degree of gelation.

Figure 5

Cross-linking probability for different molecular weight, system size, and chemical potential of dynamic bonds. The collapse over N and V indicates that the binding equilibrium is local. The dashed line shows the best fit, p(μ) = 9.57–1.11μ + 0.03μ2.

As a preliminary step toward quantifying the degree of gelation, we assume that the effect of loops is negligible, and that the system reaches the sol–gel transition when each chain has, on average, one outgoing cross-link. This gives the mean-field estimate to the gel point, . The true gel point pc will be estimated below. It is generally higher than p0 because the formation of loops or other defects reduces the effective cross-linking, requiring a higher value of p to produce a percolating network. Without prior knowledge of pc, we use the mean-field value p0 to define the degree of gelationwhich quantifies the distance from the mean-field gel point.

Cluster Number Distribution

We ran the cluster number distribution simulations for molecular weights N = 12, 25, 50, and system size V/V0 = 1, 8, 32, 64. For each N and V, we collected 11 sets of simulations corresponding to ϵ values between 0 and 1.0 with an interval of 0.1. The sets of μ values are different for different N, because p0 varies with N. Choosing nearly identical ϵ values ensures that the systems have a similar degree of gelation.

For a given p or ϵ value, the cluster number distribution function n(s) is calculated as the number of clusters containing s precursor chains divided by the total number of chains. The value of s ranges from 1 to the size of the largest cluster observed. Clusters that traverse the simulation box are unfolded according to the periodic boundary condition. The infinite clusters found past the gel point are those percolating the box along any dimension, whose s value is constrained by the finite simulation box. To improve statistics, for the same N and V, we have averaged n(s) over periodically sampled configurations.

Figure shows the distribution n(s) on a log–log scale for different N and ϵ sampled from the largest system size (the corresponding results in the smaller systems are included as Figures S2–S4 in the Supporting Information and compared in Figure S5). A remarkable agreement with the trend expected from the scaling ansatz, eq , is found. For ϵ ≤ 0.4, the cluster number decreases with s according to a power law, until reaching a cutoff value s*, followed by a rapidly decaying tail. The cutoff s* is the characteristic cluster number for the typical cluster formed at a given ϵ.[19] The physical dimension or the radius of gyration for clusters containing s* precursor chains is the correlation length of the system.[18,19] For p < pc, the value of s* increases with ϵ, suggesting that a greater number of large clusters are formed as the degree of gelation increases.

Figure 6

Cluster number distributions n(s) for different simulations with ϵ ∈ [0, 1] are plotted here on a log–log scale. The characteristic power law scaling of n(s) is cut off at a characteristic cluster size s*, above which the cluster numbers decay to zero. This is most apparent for ϵ near 0 and 1. For ϵ ≈ 0.5, the cutoff size appears to diverge, and the distribution follows an almost perfect power law, indicating that the system is near pc. For ϵ > 0.5, the presence of a percolating cluster is indicated by a sharp peak at large s.

At ϵ = 0.5, the cutoff is barely visible for all three N values. The data nearly follow a single power-law decay, indicating the absence of a characteristic length scale. The system is extremely close to the gel point, and we may deduce that ϵc is about 0.5. More detailed discussion of the slope and the location of the gel point will be provided below.

For ϵ ≥ 0.6, the cutoff s* reappears, and its value decreases with ϵ, fully compatible with the expectation that the correlation length decreases further away from the gel point. Moreover, a well-defined peak separated from the power-law/cutoff regime emerges in the region of large s values. These peaks represent the percolating clusters or networks formed past the gel point. The location of these peaks increases with ϵ, suggesting that more precursor chains have been connected to the networks. For ϵ = 1.0, the number of precursor chains in these networks become comparable to the total number of precursor chains in the system (about 3000, 6000, 12000 for N = 50, 25, 12), suggesting nearly complete gelation.

The above discussion implies that the actual gel points for all three molecular weights are around ϵ = 0.5. The average cluster size will be analyzed in the following section to estimate the gel point pc. The data from the smaller boxes (Figures S1–S5 in the Supporting Information) are entirely consistent with those from V = 64V0, except that the cluster numbers for the largest s values exhibit a hump due to finite size effects for the data near the gel point. For the systems with degrees of gelation further away from ϵ = 0.5, the cluster number distributions in Figure are not contaminated by the box size.

Gel Point

We estimate the gel point by examining the variation of the second moment of the cluster number distribution, defined as . By our convention, the denominator gives the fraction of sol chains, which is identity below the gel point, and which decreases with ϵ above the gel point. The fraction of gel chains is complementary, given by . The two summations in the expression of M2 include contributions from all the clusters below the gel point, but excludes those from the percolating clusters above the gel point. Because the distribution peaks (Figures S2–S5 in the Supporting Information) representing the percolating clusters are well-separated from the finite clusters, for all three molecular weights with V = 64V0, we choose convenient thresholds, s = 2000, 1000, and 500 for N = 12, 25, and 50, respectively. The percolating clusters are those above the threshold values. The data around these thresholds are on the order of 10–6 and are statistically insignificant. The data in Figure S6 show that varying these threshold values does not affect the scaling of the second moment.

The variation of M2 with ϵ is shown in the inset of Figure a for all three N values. The data around ϵ = 0.5, a crude estimate to the gel point, are excluded because the division between finite and percolating clusters is blurred. The value of M2 increases from both sides as the gelation transition is approached, and appears to diverge. This divergence is expected in the neighborhood of the gel point, and a power-law behavior M2 ∝ |p – pc|–γ holds on both sides with a common exponent γ.[18] For convenience, in place of pc, we shall use ϵc to denote the gel point; the scaling behaviors are not affected because p and ϵ are linearly related. Following literature,[19,39] we adjust the value of ϵc so that the identical power-law behaviors are found both below and above the gel point. Choosing ϵc = 0.53, 0.66, and 0.69 for N = 12, 25, and 50 results in the collapse of data over all N values shown in Figure . The common slopes give an estimate to the exponent γ = 1.62, which compares favorably to the tabulated value 1.80 for critical percolation.[19] Furthermore, the ratio of the amplitudes to the two power laws below and above the gel point is 1.8, consistent with recent results.[40] Finally, the apparent increment of ϵc with N does not imply that greater degree of cross-linking is needed. Once the normalization factor p0 in eq is restored, the probability pc decreases with N (Figure b).

Figure 7

Scaling of the second moment and the gel fraction. (a) Variation of the second moment M2 with ϵ, below and above the gel point. Inset: same data on a linear scale. (b) Variation of the gel fraction Pgel with ϵ. Inset: same data on a linear scale.

Figure 11

Cause of the peak in the cutoff function. (a) Schematic for the evolution of cluster number distribution with ϵ. The cutoff location and the amplitude are both constrained by mass conservation. (b) Variation of the gel point with N. The mean-field estimate to gel point: p0; the gel point estimated from the scaling of the second moment M2: pc. The values of p0, pc and their difference decrease with N.

The gel fraction follows a similar scaling pattern near the gel point. The inset of Figure b shows Pgel against the difference ϵ – ϵc for the three N values, using the ϵc values estimated from M2 (results from alternative threshold values, shown as Figure S7 in the Supporting Information, are entirely consistent). The sharp transition is rounded near the gel point due to the finite size effects, which are more severe for the cases with N = 25 and 50. Percolation models expect that near the gel point.[18,19]Figure shows that the data are consistent with this scaling, giving an exponent β = 0.6 for N = 25 and 50, and β = 0.44 for N = 12. Comparing these exponents to the tabulated value β = 0.41,[19] we conclude that our systems are more aligned with critical percolation. The greater discrepancy in systems with N = 25 and 50 is a result of more severe finite size effects.

The scaling of the characteristic cluster number s* can be examined as follows. Below the gel point, the value of s* may be obtained by scaling the cluster size s with s*, whose value is chosen to yield a collapse of the function n(s)sτ. The quality of such a collapse is shown in the inset of Figure . The mean-field exponent τ = 2.5 was used. However, given the small difference between this value and the critical exponent τ = 2.18, the same scaling behavior to s* was obtained. The variation of s* with ϵ from the three N values fits to the same power-law relation, , with 1/σ = 2.8 and σ = 0.36, which compares favorably to the critical exponent σ = 0.45.[19]

Figure 8

Characteristic cluster number for different ϵ values. Common scaling behaviors are obtained for different N values.

In short, the above results demonstrate that the same set of ϵc values gives consistent scaling behaviors for the second moment, the gel fraction, and the characteristic cluster number for three different precursor chain lengths, and that the exponents are more compatible with critical than mean-field scaling.

Cutoff Function

The scaling ansatz states that the power-law behavior in the cluster number is followed up to s = s*, beyond which a rapid decay is captured by a cutoff function f(s/s*). A general argument based on mass conservation[19] suggests that the cutoff function has a peak below the gel point and decays above the gel point. In this section, we try to extract the shape of the cutoff function.

Having estimated the gel point, we examine the departure from the power-law scaling in three regimes separately: well below the gel point (Figure a), well above the gel point (Figure b), and near the gel point (Figure c). Although the results from N = 12 are discussed because they are less susceptible to the finite size effects, the trends from N = 25 and N = 50 are entirely consistent.

Figure 9

Cutoff function for N = 12 at a range of ϵ values: (a) ϵ = 0 (circle) and 0.1 (triangle); (b) ϵ = 0.9 (circle) and 1.0 (triangle); (c) ϵ = 0.4 (circle) and 0.5 (triangle). (d) Cutoff function combining data below (z < 0) and above (z > 0) the gel point.

Figure a shows the variation of ns5/2 with s for ϵ = 0 and 0.1, where 5/2 is the mean-field value for the exponent τ. The peak expected for the cutoff function below the gel point is clearly seen. The peak position is around s = 6 for ϵ = 0 and s = 12 for ϵ = 0.1. The shift in the peak position is expected since the characteristic cluster number, i.e., the onset of the departure from the power-law scaling, increases with increasing ϵ. The inset of Figure a shows the log–linear plot for the same data, indicating that ns5/2 behaves as an exponential decay with s. The difference in the decay rate again reflects the difference in the characteristic cluster number s*.

Figure b shows ns5/2 for ϵ = 0.9 and ϵ = 1, two values well above the gel point. A monotonic decay is observed, with greater decay rates for the larger ϵ value. The inset shows again that the product ns5/2 in this regime fits to an exponential decay.

Figure c shows ns2.18 for ϵ = 0.4 and ϵ = 0.5, two values close to the gel point ϵc = 0.53. Note that, unlike the above two cases, the critical exponent τ = 2.18 is used for these two sets of data. The plateauing behavior around 0.2 implies that, for these two ϵ values, the cluster number n follows the critical scaling n ≃ 0.2s–2.18. Therefore, the values ϵ = 0.4 and ϵ = 0.5 fall within the critical regime.

One implication of the above analyses is that the data for 0.2 ≤ ϵ ≤ 0.3 and 0.6 ≤ ϵ ≤ 0.8 will be in the transition regime from the mean-field to critical scalings. The value of n(s) scales as s–5/2 for small s, scales as s–2.18 for intermediate s, and decays rapidly for large s. The transition between the first two regimes is at the Ginzburg point,[18]. The transition between the last two regimes is at s* ∝|ϵ – ϵc|–1/σ, with σ = 0.45.[19] Because the prefactor for these two transition points is unknown, and because the transition is broad (Figure ), it is difficult to identity these two transition points unambiguously. We therefore only use the data from the mean-field regimes to construct the crossover function.

For N = 12, we take the data in Figure a,b and define the cutoff function f(z) as follows.[19] First, we introduce the ratio , which is negative below the gel point and positive above the gel point. The value s* is chosen such that it leads to the collapse of the data in Figure a,b. The choice of exponent σ only affects the continuity at z = 0 but does not change the shape of the cutoff function. For convenience, we choose the mean-field value σ = 0.5 (see discussion below on how this choice affects the z-dependence). Second, the function value f(z) is given by the product ns5/2 for the corresponding s value. The cutoff function as defined is plotted against z in Figure d, which has exactly the same shape found from simulations for lattice percolation.[19]

The literature has postulated a Gaussian form for the cutoff function.[41] The attempted Gaussian fit in Figure a suggests that this form, although capturing the primary feature around the peak, fails to capture the slower decay further away from the peak. Finally, we repeated the above analyses for N = 25 and 50 and found that, surprisingly, the same shape of the cutoff function is obtained, as shown in Figure b. The inset further indicates that, in the regime above the gel point, the cutoff function decays exponentially with s/s*. The decay rate with z would depend on which exponent σ is used to convert s/s* to z.

Figure 10

Shape of the cutoff function. (a) Comparison of data for N = 12 and a Gaussian fit. (b) Collapse of data from different N values. N = 12: ϵ = 0, 0.1, 0.9, 1.0; N = 25: ϵ = 0, 0.1; N = 50: ϵ = 0, 0.1, 0.2.

Peak in Cutoff Function

The salient feature of the cutoff function is the presence of a peak below the gel point, which is essential for resolving the apparent inconsistency among several scaling relations,[19] and can be attributed to mass conservation. The origin of the peak has been addressed[19] and was shown to derive from the algebraic property of the condition . Here we show, alternatively, that the two-parameter scaling ansatz eq inevitably implies a single peak in the cutoff function. Since the argument to the cutoff function can be explicitly written as z = sσ(p – pc)/pc, one way to reveal the shape of f(z) is to vary p while holding s fixed. For p < pc, increasing the value of p effectively moves the value of z toward the origin.

Figure a illustrates the evolution of n(s) for a range of p values below pc on a log–log scale. The limiting behavior at pc is shown as a dashed line, and is can be referred to as nc(s). Since f(z) is proportional to the ratio n(s)/nc(s), we may focus on the shape of n(s) for arbitrary p. For all three p values in the schematic, the cluster numbers in the small s regime follow a similar slope as nc(s), whereas they all decay to zero beyond certain finite value s*. The value of s* increases with p, indicating the broadening of the power-law scaling.

The broadening requires that the amplitude of n(s) curves decreases with increasing p. Because of mass conservation , when s* increases, the value of n(s) in the small s regime has to decrease monotonically with p. Then following the vertical line representing a constant s and tracing the intersection with the n(s) curves, we find that, the value of n(s) at p2 is greater than that at both p1 and p3. Therefore, n(s) exhibits a peak in this range of p values, which translates to the peak in the cutoff function f(z). We emphasize that this is a combined result of the two parameter scaling for n(s) and mass conservation. The algebraic argument in ref.[19] shows convincingly that the nonmonotonic behavior in f(z) is unavoidable, while not excluding the possibility of multiple peaks. Here we note that, by monitoring the crossover between the power-law scaling and the cutoff regimes closely, the peak is unique: for either p < p1 or p > p3 below the gel point, the value of f(z) decreases.

The above analysis does not rely on the specific value of the exponents τ and σ, thus the conclusion that a peak is present holds irrespective of whether the cluster number is dominated by mean-field or critical scaling. We obtained the cutoff function from simulation data falling inside the mean-field regime. The data sufficiently near the gel point are needed to access the shape of the cutoff function in the critical regime, which is challenging because of severe finite size effects. The collapsed data of Figure b from different N values is encouraging. However, to fully demonstrate the independence of the cutoff function on molecular weight (in the mean-field regime), simulations of higher molecular weight may be needed. Increasing the molecular weight widens the mean-field scaling ranges in both p and s. One piece of evidence for this scaling is shown in Figure b for the variation of the gel point with N. The departure of the gel point from the mean-field estimate decreases with N, consistent with the expectation that the mean-field (Flory–Stockmayer) theory becomes asymptotically accurate as molecular weight increases, due to the enhanced overlap of clusters for large N.[40]

Summary

This work focuses on the distribution of dynamic clusters near the gel point. Our main result is the cutoff function, obtained from data for three N values, falling within the mean-field scaling regimes, both below and above the gel point. The gel points were estimated by studying the scaling of the second moment, the gel fraction, and the characteristic cluster number. The shape of the cutoff function is similar to that obtained from the simulation of the site percolation problem on a 2D lattice. A Gaussian fit to the cutoff function works well below the gel point but overestimates the decay rate above the gel point. Although the precise dependence of f(z) on z would depend on which σ value is used to relate z to s/s*, our data suggest that the dependence on s/s* is exponential on both sides of the gel point. We stress, however, that the cutoff function was constructed from data following mean-field scaling. Systems extremely close to the gel point have diverging correlation lengths and suffer more severe finite size effects. Although not able to identify the shape of the cutoff function when the cluster number distributions are dominated by the critical scaling, we speculate that the shape of the cutoff function is analogous. In future work, we will present the results on the cluster size, degree of overlap, loops or defect statistics, and lifetime of dynamic clusters.