Generalized Deam–Edwards approach to the statistical mechanics of randomly crosslinked systems

The crosslinking process that we shall consider is an instantaneous crosslinking process, as already mentioned previously. Many springs (crosslinkers) are introduced into the liquid system (in its preparation state), completely randomly and independently of one another. We further assume that the extension of every spring d is distributed according to a Gaussian equilibrium Gibbs–Boltzmann factor with characteristic lengthscale b (cf equation (4)):

$$\begin{equation} P( \boldsymbol{d} ) \propto \mathrm{e}^{ -| \boldsymbol{d} |^2/{2 b^2} }. \end{equation} \tag{ 8 }$$

In addition, we assume that each functional site c⁰_i has an 'effective region' of radius , centered around c⁰_i. If an end of a spring is introduced at some position x that lies within the 'effective region' around c⁰_i, it is successfully linked to this functional site. We assume that is smaller than half the minimal distance between any two functional sites, so that any end of a spring cannot be linked to two sites simultaneously. This assumption does not limit the generality of our approach (figure 2).

Zoom In Zoom Out Reset image size

Next, let p_k(c⁰_i,c⁰_j) be the probability that a given pair of sites {c⁰_i,c⁰_j} is linked by some integer number k of springs. If we assume that the crosslinkers are introduced into the system at random, independently and homogeneously, then k would be Poisson-distributed

$$\begin{equation} p_{k}({\boldsymbol{c}}^0_i,{\boldsymbol{c}}^0_j) = C \, \frac{\hat{\mu}^k}{k!}\, \mathrm{e}^{- k\, \left| {\boldsymbol{c}}^0_i - {\boldsymbol{c}}^0_j \right| ^2/2 b^2}, \end{equation} \tag{ 9 }$$

where we have introduced the dimensionless control parameter $\hat {\mu }$ for the crosslinking density⁶, and the constant C is determined via overall normalization, see below. Note that p_k(c⁰_i,c⁰_j) contains a factor proportional to the k^th power of Gibbs–Boltzmann factor (8), which characterizes the equilibrium distribution of each spring. It is convenient to define a (three-dimensional) soft delta function δ_b(x) in terms of a Gaussian spatial profile having variance parameter b²

$$\begin{equation} \delta_b({\boldsymbol{x}}) \equiv ( \sqrt{2\pi}/b )^{3} \,\mathrm{e}^{- {\boldsymbol{x}}^2/2 b^2}. \end{equation} \tag{ 10 }$$

The function δ_b(x) tends to the Dirac delta function as b → 0. The linking probability (9) can then also be expressed in the following concise form:

$$\begin{equation} p_{k}({\boldsymbol{c}}^0_i,{\boldsymbol{c}}^0_j) = \mathrm{e}^{- \tilde{\mu} \, \delta_b({\boldsymbol{c}}^0_i - {\boldsymbol{c}}^0_j)} \frac{1}{k!} \left(\, \tilde{\mu} \, \delta_b({\boldsymbol{c}}^0_i - {\boldsymbol{c}}^0_j) \,\right)^k, \end{equation} \tag{ 11 }$$

where we have exchanged the crosslink density control parameter $\hat {\mu }$ for the parameter $\tilde {\mu }$, defined via

$$\begin{equation} \tilde{\mu} \equiv (b/ \sqrt{2\pi} )^{3} \, {\hat{\mu}}. \end{equation} \tag{ 12 }$$

It is easy to check that equation (11) is properly normalized: $\sum _{k = 0}^{\infty } p_k({\boldsymbol {c}}^0_i,{\boldsymbol {c}}^0_j) =1$.

We assume that the liquid, prior to crosslinking, is in thermal equilibrium with respect to Hamiltonian H⁰_liq(c⁰), where c⁰ ≡ {c⁰_i} denotes the locations of all N functional sites in the preparation state⁷. The joint probability density function for all N sites (prior to crosslinking) to be at the locations {c⁰₁,...,c⁰_N} is then given by

$$\begin{equation} P({\boldsymbol{c}}^0) = \frac{1}{Z^0_{\rm{liq}} }\, \mathrm{e}^{- H^0_{\rm{liq}}({\boldsymbol{c}}^0)}, \end{equation} \tag{ 13 }$$

where the normalizing partition function is given by

$$\begin{equation} Z^0_{\rm{liq}} = \int D {\boldsymbol{c}}^0 \, \mathrm{e}^{- H^0_{\rm{liq}}({\boldsymbol{c}}^0)}. \end{equation} \tag{ 14 }$$

Now recall that the network structure is characterized by a set of N(N + 1)/2 stochastic variables⁸, or link numbers, {k_(i,j)}, each of them characterizing the number of springs linking a particular pair (i,j). Given a liquid configuration c⁰ prior to crosslinking, these link numbers are fully uncorrelated, and governed by the following conditional probability distribution:

$$\begin{eqnarray} P({\chi}|{\boldsymbol{c}}^0) &=& \prod_{(i,j)}p_{k_{(i,j)}}({\boldsymbol{c}}^0_i, {\boldsymbol{c}}^0_j)\nonumber\\ &=& \prod_{(i,j)}\mathrm{e}^{- \tilde{\mu} \, \delta_b({\boldsymbol{c}}^0_i - {\boldsymbol{c}}^0_j)} \frac{1}{k_{(i,j)}!} \left( \tilde{\mu} \, \delta_b({\boldsymbol{c}}^0_i - {\boldsymbol{c}}^0_j) \right)^{k_{(i,j)}}, \end{eqnarray} \tag{ 15 }$$

where we have used the linking probability given by equation (11). Making the further definition

$$\begin{equation} H_{\rm{norm}} ({\boldsymbol{c}}) \equiv - \tilde{\mu} \, \sum_{i, j } \delta_b({{\boldsymbol{c}}}_i - {{\boldsymbol{c}}}_j), \end{equation} \tag{ 16 }$$

we can write this conditional probability as

$$\begin{eqnarray} P({\chi}|{\boldsymbol{c}}^0) &=& \mathrm{e}^{H_{\rm{norm}}({\boldsymbol{c}}^0)} \prod_{(i,j)} \frac{1}{k_{(i,j)}!} \left( \tilde{\mu} \, \delta_b({\boldsymbol{c}}^0_i - {\boldsymbol{c}}^0_j) \right)^{k_{(i,j)}}\nonumber\\ &= & \mathrm{e}^{H_{\rm{norm}}({\boldsymbol{c}}^0)} \prod_{(i,j)} \frac{\hat{\mu}^{k_{(i,j)}}}{k_{(i,j)}!} \mathrm{e}^{- k_{(i,j)}\,|{\boldsymbol{c}}^0_i - {\boldsymbol{c}}^0_j|^2/2 b^2}\nonumber\\ &=& \mathrm{e}^{H_{\rm{norm}}({\boldsymbol{c}}^0)-\Delta H_{\chi}({\boldsymbol{c}}^0)} \prod_{(i,j)} \frac{\hat{\mu}^{k_{(i,j)}}}{k_{(i,j)}!}, \end{eqnarray} \tag{ 17 }$$

where we have used equations (10), (12) and (4) in the above three equalities, respectively. If we were to specify the network structure using the set χ' = {(i_e,j_e),e = 1,...,M} instead of the set χ, the corresponding conditional probability would be

$$\begin{equation} P(\chi'|{\boldsymbol{c}}^0) = \mathrm{e}^{H_{\rm{norm}}({\boldsymbol{c}}^0)} \frac{1}{M!} \prod_{e = 1}^M \tilde{\mu} \, \delta_b({\boldsymbol{c}}^0_{i_e} - {\boldsymbol{c}}^0_{j_e}), \end{equation} \tag{ 18 }$$

where we have divided the probability by the proper degeneracy factor (3), so as to cancel the over-counting.

The probability of obtaining the network structure χ can be arrived at by using the law of total probability as

$$\begin{eqnarray} P_{\chi} &=&\sum_{{\boldsymbol{c}}^0} P(\chi | {\boldsymbol{c}}^0) P({\boldsymbol{c}}^0)\nonumber\\ &=& \frac{1}{Z^0_{\rm{liq}}} \int D{\boldsymbol{c}}^0 \, \mathrm{e}^{- H_{\rm{liq}}^0({\boldsymbol{c}}^0) + H_{\rm{norm}}({\boldsymbol{c}}^0) -\Delta H_{\chi}({\boldsymbol{c}}^0)}\, \prod_{(i,j)} \frac{\hat{\mu}^{k_{(i,j)}}}{k_{(i,j)}!}\,. \end{eqnarray} \tag{ 19 }$$

Furthermore, by making the definition

$$\begin{equation} \tilde{Z}_{0,\chi} \equiv \int D {\boldsymbol{c}}^0 \, \mathrm{e}^{-H_{\rm{liq}}^0({\boldsymbol{c}}^0) + H_{\rm{norm}}({\boldsymbol{c}}^0) -\Delta H_{\chi}({\boldsymbol{c}}^0)}\,, \end{equation} \tag{ 20 }$$

we obtain

$$\begin{equation} P_{\chi} = \frac{\tilde{Z}_{0,\chi}}{Z^0_{\rm{liq}}} \prod_{(i,j)} \frac{\hat{\mu}^{k_{(i,j)}}}{k_{(i,j)}!}. \end{equation} \tag{ 21 }$$

Hence we see that, apart from a dimensionless factor, P_χ is the ratio of two partition functions. Even though the variables c⁰_i (associated with preparation state) appear only as dummy variables, the integrations over them in equation (20) cannot be carried out explicitly. Below, when we calculate the disorder average of free energy, we shall find that these variables c⁰_i interact with the corresponding variables c_i in the measurement ensemble. It is this interaction that generates the ubiquitous memory effects in randomly crosslinked materials.

The conditional probability that there are precisely M linking springs in the system can be obtained by summing equation (18) over all possible values taken by the indices (i_e,j_e):

$$\begin{eqnarray} P(M|{\boldsymbol{c}}^0) &=&\sum_{i_1, j_1 = 1}^N \sum_{i_2, j_1 = 2}^N\cdots\sum_{i_M, j_M = 1}^NP(\chi' | {\boldsymbol{c}}^0)\nonumber\\ &=& \mathrm{e}^{H_{\rm{norm}}({\boldsymbol{c}}^0)}\frac{1}{M!}\left(- H_{\rm{norm}}({\boldsymbol{c}}^0)\right)^{M}. \end{eqnarray} \tag{ 22 }$$

It is elementary to check that the conditional probability P(M|c₀) is properly normalized for any fixed c₀:

$$\begin{equation} \sum_{M = 0}^{\infty} P(M | {\boldsymbol{c}}^0) = 1. \end{equation} \tag{ 23 }$$

To check the normalization of the probability P_χ, we sum equation (19) over all possible network structures. By equation (2), this amounts to summing over all possible linking numbers {k_(i,j) }. Further using equations (4) we arrive at

$$\begin{eqnarray} \sum_{\chi} P_{\chi} &=&\frac{1}{Z^0_{\rm{liq}}} \int D{\boldsymbol{c}}^0 \,\mathrm{e}^{- H_{\rm{liq}}^0({\boldsymbol{c}}^0) + H_{\rm{norm}}({\boldsymbol{c}}^0)}\sum_{\chi} \mathrm{e}^{-\Delta H_{\chi}({\boldsymbol{c}}^0)} \,\prod_{(i,j)} \frac{\hat{\mu}^{k_{(i,j)}}}{k_{(i,j)}!}\nonumber\\ &= &\frac{1}{Z^0_{\rm{liq}}} \int D{\boldsymbol{c}}^0 \, \mathrm{e}^{- H_{\rm{liq}}^0({\boldsymbol{c}}^0) + H_{\rm{norm}}({\boldsymbol{c}}^0)} \prod_{(i,j)} \sum_{k = 1}^{\infty}\frac{\hat{\mu}^k}{k!} \mathrm{e}^{ - k \left| {\boldsymbol{c}}_i - {\boldsymbol{c}}_j \right|^2/2 b^2}. \end{eqnarray} \tag{ 24 }$$

The summation and product can readily be calculated, leading to the factor $\exp - H_X({\boldsymbol {c}}^0)$, which precisely cancels the corresponding exponential factor. Hence, as expected, we have the normalization condition

$$\begin{equation} \sum_{\chi} P_{\chi} = 1. \end{equation} \tag{ 25 }$$

We are now ready to compute the disorder-averaged free energy, equation (7). For this purpose, we use the following identity, which is commonly called the replica trick⁹:

$$\begin{equation} \overline {\ln X} = \lim_{n \rightarrow 0} \frac{1}{n} \ln \overline{X^n}. \end{equation} \tag{ 26 }$$

We use this to rewrite the average free energy as

$$\begin{eqnarray} \overline{ F} &=& - T \lim_{n \rightarrow 0} \frac{1}{n} \ln \overline{Z_{\chi}^n}= - T \lim_{n \rightarrow 0} \frac{1}{n} \ln\sum_{\chi} P_{\chi} Z_{\chi}^n \nonumber\\ &=& - T \lim_{n \rightarrow 0} \frac{1}{n} \ln\sum_{\chi} \left( \prod_{(i,j)} \frac{\hat{\mu}^{k_{(i,j)}}}{k_{(i,j)}!} \right) \frac{1}{Z_{\rm{liq}}^0}\,\, \tilde{Z}_{0,\chi}\,(Z_{\chi})^n. \end{eqnarray} \tag{ 27 }$$

We calculate the right-hand side of equation (27) for positive integers n, and subsequently analytically continue n to real values. More specifically, we calculate Zⁿ_χ by replicating the variables {c_i} a number n times, to obtain an n-tuple (c¹_i,c²_i,...,cⁿ_i). Together with the corresponding variable c⁰_i of the preparation ensemble, they form a (1 + n)-tuple $\hat {{\boldsymbol {c}}}_i \equiv ({\boldsymbol {c}}^{0}_i, {\boldsymbol {c}}^1_i,\ldots, {\boldsymbol {c}}^n_i )$. The shorthand $\hat {{\boldsymbol {c}}}$ has been extensively used for this (1 + n)-tuple in the vulcanization theory literature. Using equations (6) and (20), the product of partition functions Z_0,χ(Z_χ)ⁿ in equation (27) can be written in terms of functional integrals (with the shorthand notation $D\hat {{\boldsymbol {c}}} \equiv \prod _{\alpha = 0}^n D {\boldsymbol {c}}^{\alpha }$) as follows:

$$\begin{equation} Z_{0,\chi}\,( Z_{\chi})^n = \int D \hat{{\boldsymbol{c}}} \, \prod_{\alpha = 0}^n \, \mathrm{e}^{- \tilde{H}^0_{\rm{liq}} ({\boldsymbol{c}}^0) - \sum_{\alpha = 1}^n H_{\rm{liq}} ({\boldsymbol{c}}^{\alpha}) - \sum_{\alpha = 0}^n \Delta H_{\chi}({\boldsymbol{c}}^{\alpha}) }. \end{equation} \tag{ 28 }$$

Next, by using the definition equation (4) we easily see that

$$\begin{eqnarray} \sum_{\alpha = 0}^n \Delta H_{\chi}({\boldsymbol{c}}^{\alpha}) &=& \sum_{\alpha = 0}^n \frac{1}{2 b^2} \sum_{(i,j)} k_{(i,j)} \left| {{\boldsymbol{c}}}_{i} - {{\boldsymbol{c}}}_{j} \right| ^2 \nonumber\\ &=&\frac{1}{2 b^2} \sum_{(i,j)} k_{(i,j)} \left| {\hat{{\boldsymbol{c}}}}_{i} - \hat{{\boldsymbol{c}}}_{j} \right|^2 \nonumber\\ &\equiv& \Delta H_{\chi} (\hat{{\boldsymbol{c}}}), \end{eqnarray} \tag{ 29 }$$

where we have used the notation $|\hat {{\boldsymbol {c}}}^{\alpha }|^2 \equiv \sum _{\alpha = 0}^n |{\boldsymbol {c}}_i^{\alpha }|^2$. We can now substitute equation (28) back into equation (27), exchange the order in which $\sum _{\chi }$ and $\int D \hat {{\boldsymbol {c}}}$ are performed in equation (27), and sum over χ for a fixed configuration of $\hat {{\boldsymbol {c}}}$. The summation over the realizations of network structure can be calculated using the following result:

$$\begin{eqnarray} &&\fl \sum_{\chi} \mathrm{e}^{- \Delta H_{\chi}(\hat{{\boldsymbol{c}}})} \left( \prod_{(i,j)} \frac{\hat{\mu}^{k_{(i,j)}}}{k_{(i,j)}!} \right)= \sum_{\{k_{(i,j)} \}} \left( \prod_{(i,j)} \frac{\hat{\mu}^{k_{(i,j)}}}{k_{(i,j)}!} \right) \mathrm{e}^{- \frac{1}{2 b^2} \sum_{(i,j)} k_{(i,j)} \left| {\hat{{\boldsymbol{c}}}}_{i} - \hat{{\boldsymbol{c}}}_{j} \right|^2} \nonumber\\ &&= \prod_{(i,j)} \sum_{k} \left( \frac{\hat{\mu}^{k}}{k!} \right) \mathrm{e}^{- \frac{1}{2 b^2} k \left| {\hat{{\boldsymbol{c}}}}_{i} - \hat{{\boldsymbol{c}}}_{j} \right|^2} \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} &=&\prod_{(i,j)} \exp \left( \hat{\mu} \, \mathrm{e}^{- \frac{1}{2 b^2} k \left| {\hat{{\boldsymbol{c}}}}_{i} - \hat{{\boldsymbol{c}}}_{j} \right|^2} \right)\nonumber\\ &=& \exp \sum_{(i,j)} {\mu} \, \delta_b(\hat{{\boldsymbol{c}}}_i - \hat{{\boldsymbol{c}}}_j), \end{eqnarray} \tag{ 31 }$$

where in the last equality, we have used the definitions

$$\begin{equation} {\mu} \equiv \hat{\mu} (\sqrt{2\pi}/b)^{3(1+n)}, \end{equation} \tag{ 32 }$$

$$\begin{equation} \delta(\hat{{\boldsymbol{x}}}) \equiv \prod_{\alpha = 0}^n \delta_b({\boldsymbol{x}}^{\alpha}) = \delta_b({\boldsymbol{x}}^0) \delta_b({\boldsymbol{x}}^1) \cdots \delta_b({\boldsymbol{x}}^n). \end{equation} \tag{ 33 }$$

As we shall eventually take the replica limit, n → 0, the difference between μ and $\tilde {\mu }$ (see equation (12)) can be ignored.

Combining all these results, we find that the disorder-averaged free energy is given as follows:

$$\begin{equation} \overline{F} = - T\, \lim_{n \rightarrow 0} \frac{1}{n} \, \ln\left( \frac{1}{Z^0_{\rm{liq}}} \, Z_{1+n}\right), \end{equation} \tag{ 34a }$$

$$\begin{equation} Z_{1+n} = \int D\hat{{\boldsymbol{c}}} \, \mathrm{e}^{-H_{1+n}(\hat{{\boldsymbol{c}}}) }, \end{equation} \tag{ 34b }$$

$$\begin{equation} H_{1+n}(\hat{{\boldsymbol{c}}}) = H_{\rm{liq}}^0({\boldsymbol{c}}^0) - H_{\rm{norm}}({\boldsymbol{c}}^0) + \sum_{\alpha = 1}^n H_{\rm{liq}} ({\boldsymbol{c}}^{\alpha}) + H_{\rm{norm}} (\hat{\boldsymbol{c}}), \end{equation} \tag{ 34c }$$

$$\begin{equation} H_{\rm{norm}} (\hat{\boldsymbol{c}}) = - {\mu} \sum_{(i, j)}^N \delta_b \left( \hat{{\boldsymbol{c}}}_i - \hat{{\boldsymbol{c}}}_j\right). \end{equation} \tag{ 34d }$$

$H_{\rm{norm}}[\hat {\boldsymbol {c}}]$

$\hat {\boldsymbol {c}} = \{ {\boldsymbol {c}}^{\alpha } \}_{n = 0}^n$

We can use Z_1+n with appropriate source terms to calculate the average of a physical quantity A(c⁰) (≡ A⁰) that only depends on the degrees of freedom c⁰ of the preparation ensemble. As replicas 1,...,n play no role in such a calculation, we can take the replica limit n → 0 before the calculation of the expectation value. This is equivalent to deleting all variables c^α for α = 1,...,n from the replicated Hamiltonian (34c). In particular, we have

$$\begin{eqnarray*} H_{\rm{norm}}(\hat{{\boldsymbol{c}}}) &\rightarrow& H_{\rm{norm}}({\boldsymbol{c}}^0), \nonumber\\ H_{1+n}(\hat{{\boldsymbol{c}}}) &\rightarrow& H_{\rm{liq}}^0({\boldsymbol{c}}^0). \end{eqnarray*}$$

We also have to delete the integrals over c^α (for α = 1,...,n) in the functional integrals (34b) and (35). Therefore, we have

$$\begin{equation} \langle A^0 \rangle = \frac{\int D \hat{{\boldsymbol{c}}} A^0\, \mathrm{e}^{- H_{1+n}}}{\int D \hat{{\boldsymbol{c}}} \, \mathrm{e}^{- H_{1+n}} } = \frac{\int D {\boldsymbol{c}}^0 \,A^0 \, \mathrm{e}^{- H^0_{\rm{liq}} ({\boldsymbol{c}}^0) } }{\int D {\boldsymbol{c}}^0 \, \mathrm{e}^{- H^0_{\rm{liq}} ({\boldsymbol{c}}^0) } } = \langle A \rangle_{\rm{liq}}^0 \equiv [ A ]_{\rm{p}}. \end{equation} \tag{ 36 }$$

As a result, 〈A⁰〉 is indeed the average in the preparation state, as expected.

In the replicated partition function, equation (34b), there is one copy of the preparation ensemble (the zeroth replica), and n copies of the measurement ensemble. In figure 3 we show a cartoon that shows schematically the relationship between various replicas. All 1 + n replicas interact with one another through the last term in equation (34c), which is linear in the crosslinking density. Note that there are n copies of the measurement ensemble (i.e. replicas 1,...,n) but only one copy of the preparation ensemble (i.e. replica 0). In general, owing to the interactions, the preparation and measurement ensembles can mutually influence each other. For example, if we tune some parameter J⁰ in fluid Hamiltonian H⁰_liq of the preparation state (see equation (34c)), there will be changes in the physical quantities in the measurement state

$$\begin{equation} \frac{\partial}{\partial J^0} \left[\overline{\langle A \rangle_{\rm{m}}}\right]_{\rm{p}} \neq 0. \end{equation} \tag{ 37 }$$

This is, of course, entirely reasonable, as the properties of polymer networks do depend on their method of preparation. Reciprocally, however, if we change a parameter J in the measurement state, might there, at least in principle, also be responses of physical quantities in the preparation state, e.g.

$$\begin{equation} \frac{\partial}{\partial J} [A ]_{\rm{p}} \neq 0 \quad (?). \end{equation} \tag{ 38 }$$

Such a response would be entirely unphysical, as it would violate the principle of causality. In the framework of vulcanization theory, the causality principle is rescued via the replica trick: as n → 0, the weight of the measurement ensemble relative to the preparation ensemble tends to zero, and therefore the former can have no quantitative influence on the latter.

Zoom In Zoom Out Reset image size

An explicit proof of this result is simple. Let the parameter J couple linearly to some quantity Φ(c) in the (1 + n)-replica Hamiltonian H_1+n; this leads to the change

$$\begin{eqnarray*} &&H_{1+n} \rightarrow H_{1+n} - J \sum_{\alpha = 1}^n \Phi({\boldsymbol{c}}^{\alpha}). \end{eqnarray*}$$

Substituting this back into the expectation value equation (36), taking the derivative with respect to J, and finally setting J = 0, we obtain

$$\begin{equation} \frac{\partial [ A ]_{\rm{p}} }{\partial J} = \sum_{\alpha = 1}^n \langle A^0 \Phi({\boldsymbol{c}}^{\alpha}) \rangle = n \langle A^0 \Phi({\boldsymbol{c}}^{1}) \rangle, \end{equation} \tag{ 39 }$$

where in the last equality we have assumed that permutation symmetry amongst replicas 1 to n remains intact. The right-hand side of equation (39) therefore vanishes in the replica limit, i.e. n → 0. Even in the presence of the spontaneous breaking of this replica permutation symmetry, the right-hand side of equation (39) would even be proportional to n, and hence would vanish in the replica limit.

In the replicated Hamiltonian (34c), the parameters of the preparation state only appear in H⁰_liq, whereas the parameters of the measurement state only appear in H_liq. These two types of parameters can be separately controlled in experiments, but they jointly determine the statistical physics of randomly crosslinked materials. Let us consider a rather special case in which H⁰_liq = H_liq, i.e. a system that is prepared and measured under identical conditions. What appears inside the (1 + n)-replica Hamiltonian (34c) is, however, H⁰_liq(c⁰) − H_norm(c⁰) instead of just H⁰_liq(c⁰):

$$\begin{equation} H_{1+n}(\hat{{\boldsymbol{c}}}) = H_{\rm{liq}}({\boldsymbol{c}}^0) - H_{\rm{norm}}({\boldsymbol{c}}^0) + \sum_{\alpha = 1}^n H_{\rm{liq}} ({\boldsymbol{c}}^{\alpha}) + H_{\rm{norm}} (\hat{\boldsymbol{c}}). \end{equation} \tag{ 40 }$$

Consequently, the replica Hamiltonian H_1+n is not invariant under permutations of the replicas that mix the zeroth and the other replicas. Hence, the average of a quantity A in the preparation state is generically different from its average in the measurement state:

$$\begin{eqnarray*} &&[ A ]_{\rm{p}} = \langle A^0 \rangle \neq \langle A^{\alpha} \rangle = \overline{ \langle A \rangle_{\rm{m}}}. \end{eqnarray*}$$

This is completely reasonable, as the random network structure does affect physical quantities in the measurement state but cannot influence the preparation state. One simple example involves the nematic correlations in liquid crystalline elastomers that have been crosslinked in isotropic state. The random polymer network tends to disorder the nematic degrees of freedom, and therefore reduces the nematic correlations.

Historically, Deam and Edwards [1] chose a particular probability distribution P_χ for the network structure given by

$$\begin{equation} P_{\chi}^{\rm{DE}} = \frac{\tilde{Z}_{0,\chi}^{\rm{DE}}}{Z^0_{\rm{liq}}} \prod_{(i,j)} \frac{\hat{\mu}^{k_{(i,j)}}}{k_{(i,j)}!}, \end{equation} \tag{ 41 }$$

where

$$\begin{equation} \tilde{Z}_{0,\chi}^{\rm{DE}} = \int D {\boldsymbol{c}}^0 \, \mathrm{e}^{-H_{\rm{liq}}^0({\boldsymbol{c}}^0) -\Delta H_{\chi}({\boldsymbol{c}}^0)}.\end{equation} \tag{ 42 }$$

The corresponding (1 + n)-replica Hamiltonian is then given by

$$\begin{eqnarray} H_{1+n}^{\rm{DE}}(\hat{{\boldsymbol{c}}}) &=&H_{\rm{liq}}({\boldsymbol{c}}^0)+ \sum_{\alpha = 1}^n H_{\rm{liq}} ({\boldsymbol{c}}^{\alpha})+ H_{\rm{norm}} (\hat{\boldsymbol{c}}), \nonumber\\ &= & \sum_{\alpha = 0}^n H_{\rm{liq}} ({\boldsymbol{c}}^{\alpha})+ H_{\rm{norm}} (\hat{\boldsymbol{c}}),\end{eqnarray} \tag{ 43 }$$

in which all 1 + n replicas appear symmetrically. The resulting 1 + n replica theory then has an enlarged S_1+n permutation symmetry, rather than the S_n permutation symmetry that is mandatory for replica theories. The Deam–Edwards choice of disorder distribution is appealing, in that it simplifies the analysis considerably. However, it also leads to the unnatural consequence that the averages of arbitrary physical quantities in the preparation state are equal to their counterparts in the measurement state. The S_1+n permutation in the Deam–Edwards theory is therefore not mandatory, as first pointed out by Broderix et al [4]. The explicit breaking of the S_1+n permutation symmetry down to S_n permutation symmetry is physically natural, and opens the way to a systematic analysis of the interrelationship between randomly crosslinked materials and their histories of formation. This is a valuable and important element of the physics contained in vulcanization theory.

We end this section with a brief discussion of the difference between the vulcanization model and certain well-known spin glass models. Consider, e.g. the Edwards–Anderson model [5]:

$$\begin{eqnarray*} &&H = - \sum_{\langle i,j\rangle} J_{ij} S_i S_j. \end{eqnarray*}$$

The quenched random variables {J_ij} are typically assumed to have independent Gaussian statistics with vanishing means, and therefore averaging over them (after application of the replica method) is elementary. The final, replicated theory reads

$$\begin{eqnarray*} &&H_{\mathrm{R}} = - \frac{1}{2} \sigma_J^2 \sum_{\langle i,j\rangle} \sum_{\alpha, \beta} S^{\alpha}_i S^{\alpha}_j S^{\beta}_i S^{\beta}_j \end{eqnarray*}$$

and contains n replicas of the annealed degrees of freedom (i.e. the spins). The variance of the quenched variables σ_J appears only as a parameter in the replicated Hamiltonian. The quenched variables {J_ij} have been integrated out, and therefore do not show up in the replica theory. By contrast, in vulcanization it is not just expedient but indeed manifestly physical to employ an equilibrium preparation state to generate the statistics for the network structure. (For vulcanization, the quenched random variables are the network structure.) Furthermore, correlations between the preparation state variables and the measurement state variables encode the permanent memory effects, and these are experimentally measurable. Finally, we note that the Edwards–Anderson model represents a considerable idealization of disordered systems in real world. In general, the physics of most disordered systems does depend on the method of preparation. Hence, the proper statistical modeling of such systems requires a separate statistical ensemble.

We can integrate out the n variables x¹,...,xⁿ in the order parameter (51), leaving only one variable x⁰, and thereby obtain a one-replica distribution

$$\begin{equation} \langle \Omega^{0}({\boldsymbol{x}}^{0}) \rangle_{\rm{eff}} \equiv \int \prod_{\alpha = 1}^n \mathrm{d} {\boldsymbol{x}}^{\alpha} \langle \Omega(\hat{{\boldsymbol{x}}}) \rangle = \left \langle \sum_{i=1}^{N} \delta_{b^{\prime}}({\boldsymbol{x}}^0 - {\boldsymbol{c}}^0_i) \right \rangle_{1+n} = \left[ \rho({\boldsymbol{x}}^0) \right]_{\rm{p}}. \end{equation} \tag{ 53 }$$

The resulting object is the average number-density of functional sites in the preparation state, which is manifestly an intensive quantity. At the coarse-grained level, it can also be understood as the average particle density of the liquid prior to crosslinking. The fluctuation δΩ⁰(x⁰) = Ω⁰(x⁰) − 〈Ω⁰(x⁰)〉_eff is also linearly related to the fluctuation of particle density in the preparation state. Similarly, we can obtain the αth replica distribution for α ≠ 0:

$$\begin{equation} \langle \Omega^{\alpha}({\boldsymbol{x}}^{\alpha}) \rangle_{\rm{eff}} \equiv \int \prod_{\beta \neq \alpha}^n \mathrm{d} {\boldsymbol{x}}^{\beta} \langle \Omega(\hat{{\boldsymbol{x}}}) \rangle = \left \langle \sum_{i=1}^{N} \delta_{b^{\prime}}({\boldsymbol{x}}^{\alpha} - {\boldsymbol{c}}^0_i) \right \rangle_{1+n} = \left[ \overline{\langle \rho({\boldsymbol{x}}^{\alpha}) \rangle_{\rm{m}} }\right]_{\rm{p}}, \end{equation} \tag{ 54 }$$

which describes average particle density in the measurement state. This relationship between the vulcanization order parameter and the particle densities in the preparation and measurement states is the primary reason that we choose the particular normalization of the order parameter field in equation (44).

For swollen gels and elastomers there can be large density fluctuations in the system (either quenched or thermal). These fluctuations, characterized by correlations of the one-replica distributions Ω^α(x^α), may strongly influence the elasticity of system. For unswollen gels and elastomers, the mean particle density at long enough lengthscales is usually uniform and its fluctuations are negligibly small. The one-replica part of free energy then contains no interesting physics, besides the stabilizing of the density fluctuations.

Finally, by integrating out n − 1 variables, we obtain the two-replica correlator

$$\begin{equation} \langle \Omega^{\alpha \beta}({\boldsymbol{x}}^{\alpha} - {\boldsymbol{x}}^{\beta}) \rangle_{\rm{eff}} = \int \prod_{\gamma \neq \alpha,\beta}^n \mathrm{d} {\boldsymbol{x}}^{\gamma} \langle \Omega(\hat{{\boldsymbol{x}}}) \rangle \end{equation} \tag{ 55 }$$

which encodes the system average of the correlation of the position fluctuations of a site in two independent measurements. Macroscopic translational symmetry dictates that this function depend only on the difference between two coordinates, x^α and x^β. Generically, it contains less information than the full order parameter $\Omega ({\hat {\boldsymbol {x}}})$.

By solving the saddle-point equation (59), subject to the constraint (57), we obtain the saddle-point value of the order parameter, which we denote by $\overline {\Omega }({\hat {\boldsymbol {x}}})$. For r > 0, the saddle-point value of $\overline {\Omega }$ is trivially zero, corresponding to the liquid state. For r < 0, the trivial saddle-point becomes unstable, and a non-trivial (i.e. replica-space dependent) saddle-point value emerges, corresponding to the amorphous solid state. It has the following form of superposed Gaussians¹⁰:

$$\begin{equation} \overline{\Omega}({\hat{\boldsymbol{x}}}) = q\, \int \mathrm{d}{\boldsymbol{z}} \int \mathrm{d} \tau \, p(\tau) \, \left({2 \pi}{\tau}\right)^{(1+n)d/2} \mathrm{e}^{- \frac{\tau}{2}\sum_{\alpha = 0}^n({\boldsymbol{x}}^{\alpha} - {\boldsymbol{z}})^2} - \frac{\,\,q }{\quad V^{n}}. \end{equation} \tag{ 60 }$$

The coefficient q is identified with the number density of the infinite cluster (i.e. the gel fraction). The dummy parameter τ is the inverse variance of the Gaussian fluctuations from their equilibrium positions of the localized particles (belonging to the infinite cluster), and is called the inverse square localization length. The fact that there is a distribution of τ signifies the heterogeneous nature of randomly crosslinked systems: the infinite cluster is more rigid in some places and looser in others. Both q and p(τ) are determined by solving the saddle-point equation. For details, see [6, 7]. It is important to note that they are independent of strain deformation and all other parameters that are tunable in the measurement state. The physics of the saddle-point order parameter can be better understood by focusing on the integrand inside equation (60). The vector z can be interpreted as the mean position of a particle in the gel fraction, and is localized around x⁰, the position of the particle at the moment of crosslinking. After crosslinking, the particle in the gel fraction (x^α,α = 1,...,n) is localized around z, its mean location in the measurement state. The uniform integral over z implies that the gel is statistically homogeneous, after averaging over the quenched disorder.

If the crosslinker density is not high enough, the system is in the sol phase in the measurement state, no infinite cluster emerges, and all particles are delocalized, so that there is no correlation between their positions in different measurements, performed in the preparation and measurement states. Hence, the associated joint pdf is a trivial constant, and the average order parameter is vanishes identically. By contrast, in the gel phase, there is an infinite cluster which constitutes a finite percentage of the overall mass. For any particle in this infinite cluster, different measurements of its position are necessarily correlated: knowledge of its location in the preparation state tells us information of its whereabout in the measurement state, and hence the average order parameter in equation (51) is non-zero. The order parameter $\Omega ({\hat {\boldsymbol {x}}})$ thus distinguishes the gel phase from the sol phase.

We can actually deduce the functional form of the saddle-point order parameter (60) using the definition of the order parameter (52) (after subtracting a trivial constant, so that the constraint (56) is satisfied) together with some simple, intuitive reasoning. We have already noted that particles belonging to the liquid fraction do not contribute to the saddle-point order parameter. For a particle i belonging to the infinite cluster, we expected that its position c_i in the measurement state is Gaussian-distributed around its mean value, which we denote by z_i, with variance 1/τ_i. The thermal averages inside equation (52) combine to give

$$\begin{equation} \langle \delta_{b^{\prime}}({\boldsymbol{x}}^1- {\boldsymbol{c}}_i) \rangle_{\rm{m}} \cdots \langle \delta_{b^{\prime}}({\boldsymbol{x}}^n- {\boldsymbol{c}}_i) \rangle_{\rm{m}} = (2 \pi \tau_i)^{nd/2}\, \exp \left\{- \frac{\tau_i}{2} \sum_{\alpha = 1}^n |{\boldsymbol{x}}^{\alpha} - {\boldsymbol{z}}_i |^2\right\}. \end{equation} \tag{ 61 }$$

The parameters z_i and τ_i of course depend on the microscopic configuration of the preparation state at the instant of crosslinking, as well as on the crosslinks that are realized. It is important to note that the mean position z_i is distinct from the position of the particle at the instant of crosslinking, which we denote by c⁰_i. Although, for a given configuration of preparation state and crosslinker realization, z_i is fully determined, it becomes a statistical variable after we include the various crosslinker realizations. In fact, we expect that z_i is Gaussian-distributed around the particle position c⁰_i in the preparation state, with the same variance 1/τ_i. Also, for different realizations of crosslinkers, the variance 1/τ_i itself should also be different, obeying some distribution p(τ). An ergodic hypothesis suggests that this distribution is reflected in the microscopic spatial heterogeneity that individual samples of elastomer exhibit. All these considerations suggest that

$$\begin{eqnarray} &&\fl \overline{ \langle \delta_{b^{\prime}}({\boldsymbol{x}}^1- {\boldsymbol{c}}_i) \rangle_{\rm{m}} \cdots \langle \delta_{b^{\prime}}({\boldsymbol{x}}^n- {\boldsymbol{c}}_i) \rangle_{\rm{m}} }= \int \mathrm{d} {\boldsymbol{z}}\, \int \mathrm{d}\tau\, (2 \pi \tau)^{(1+n)\,d/2}\,p(\tau)\nonumber\\ &&\times \exp \left\{ - \frac{\tau}{2} \sum_{\alpha = 1}^n | {\boldsymbol{x}}^{\alpha} - {\boldsymbol{z}}|^2 - \frac{\tau}{2} |{\boldsymbol{c}}_i^0- {\boldsymbol{z}}|^2 \right\}. \end{eqnarray} \tag{ 62 }$$

Finally, we note that in the preparation state the system is a homogeneous liquid and thus translation invariant. Hence, the position of a given particle c₀ is uniformly distributed over space, and therefore for an arbitrary function f(c⁰_i) we have

$$\begin{equation} \left[ f({\boldsymbol{c}}_i^0) \right]_{\rm{p}} = \int \frac{\mathrm{d}{\boldsymbol{c}}}{V}\,f({\boldsymbol{c}}). \end{equation} \tag{ 63 }$$

We now substitute equation (62) back into equation (52) and carry out the average over the preparation state and also sum over all of the particles. Thus, we finally arrive at the precise form of the saddle-point order parameter given in equation (60).

The 1 + n replicated Hamiltonian H_1+n, equation (34c), possesses the symmetry of independent translations of the replicas, i.e. it is invariant under the α-dependent translations:

$$\begin{eqnarray*} &&{\boldsymbol{c}}_i^{\alpha} \rightarrow {\boldsymbol{c}}_i^{\alpha} + {\bf u}^{\alpha}, \quad i = 1, \ldots, N. \end{eqnarray*}$$

This symmetry is also possessed by the average order parameter $\overline {\Omega }({\hat {\boldsymbol {x}}})$ in the liquid phase (which has the value zero). In the gel phase, however, this translational symmetry is explicitly broken by the saddle-point order parameter, equation (60). Translations of all particles in one replica leads to a distinct but energetically degenerate order parameter

$$\begin{eqnarray*} &&\overline{\Omega}({\boldsymbol{x}}^0, \ldots, {\boldsymbol{x}}^{\alpha} + {\boldsymbol{u}}, \ldots {\boldsymbol{x}}^n) \neq \overline{\Omega}({\boldsymbol{x}}^0, \ldots, {\boldsymbol{x}}^{\alpha}, \ldots {\boldsymbol{x}}^n). \end{eqnarray*}$$

This reduction in symmetry is a consequence of the localization of the particles associated with the emergence of an infinite cluster. Nevertheless, translations of all replicas by a common vector u remains a symmetry for equation (60):

$$\begin{eqnarray*} &&\overline{\Omega}({\boldsymbol{x}}^0,{\boldsymbol{x}}^1, \ldots, {\boldsymbol{x}}^n) = \overline{\Omega}({\boldsymbol{x}}^0 + \boldsymbol{u},{\boldsymbol{x}}^1 + \boldsymbol{u}, \ldots, {\boldsymbol{c}}^n + \boldsymbol{u}). \end{eqnarray*}$$

This macroscopic translational invariance [6] reflects the fact that gels and elastomers are statistically homogeneous. This symmetry can be broken by translational ordering, such as smectic ordering in liquid crystalline elastomers.

Because rubbery materials can typically sustain large deformations, their elasticity theory is necessarily nonlinear. This is already apparent in the classical theory of rubber elasticity, which gives the elastic free energy

$$\begin{equation} H(\boldsymbol{\Lambda}) = \frac{\mu}{2} {\rm{Tr}\,} \boldsymbol{\Lambda}^{\rm{T}} \boldsymbol{\Lambda} \end{equation} \tag{ 64 }$$

in terms of the (spatially uniform) deformation gradient matrix Λ. The matrix Λ relates the undeformed positions of particles belonging to the infinite cluster to their deformed positions z' ≡ Λ·z. As most rubbery materials are nearly incompressible, det Λ can be taken to be unity. This is largely responsible for the nonlinear nature of rubber elasticity. In the framework of statistical physics, the positions of particles fluctuate, and this deformation relation must be understood in the sense of relocation of the mean positions. It is reasonable to suppose that the saddle-point order parameter of the deformed state is given by

$$\begin{equation} \overline{\Omega}_{\boldsymbol{\Lambda}}({\hat{\boldsymbol{x}}}) = q\, \int_z \int_{\tau} \mathrm{e}^{- \frac{\tau}{2} ({\boldsymbol{x}}^0 - {\boldsymbol{z}})^2 - \frac{\tau}{2}\sum_{\alpha = 1}^n |{\boldsymbol{x}}^{\alpha} - \boldsymbol{\Lambda} {\boldsymbol{z}}|^2} - \frac{\,\,q}{\quad V^{n}}. \end{equation} \tag{ 65 }$$

Note that the mean position z is affinely deformed in the measurement state (replicas 1,...,n), but not in the preparation state (replica 0), and that q and p(τ) are not deformed (i.e. do not depend on Λ). Deformation of zeroth replica can always be 'gauged away' via a coordinate transform of the dummy integration variable z, and hence has no physical significance. It is the relative deformation between the preparation and measurement ensembles that has physical significance.

While it is rather obvious that the gel fraction q should not change with deformation, it is a priori not clear why the distribution of localization length p(τ) should stay the same in the course of deformation. Nevertheless, it has been shown explicitly [7] that equation (65) indeed satisfies the saddle-point equation (59) and the boundary conditions for a solid having a uniform deformation Λ. Furthermore, by inserting equation (65) into the Landau free energy (58)¹¹ it is straightforward to obtain the elastic free energy of the deformed system, as given in equation (64), i.e. the classical theory of rubber elasticity [8]. Moreover, one finds that μ∝|r|³ for the shear modulus¹². It is rather interesting to see that classical rubber elasticity theory emerges at the saddle-point level of vulcanization theory.

In nematic elastomers, the polymer chains carry liquid crystalline mesogens that are prone to undergoing nematic ordering. The mesogens can either be parts of the polymer backbone (i.e. main-chain nematic polymers) or attached to the backbone side-on (i.e. side-chain nematic polymer). A Landau theory of nematic elastomers has be derived via vulcanization theory elsewhere [7]; here we simply quote the results.

Let the (spatially uniform) nematic order parameter be Q⁰ in the preparation state and Q in the measurement state. The free energy (58) needs to be modified to include possible couplings between the nematic and vulcanization order parameters. To lowest order in these order parameters, the couplings are

$$\begin{equation} \int \mathrm{d} {\hat{\boldsymbol{x}}} \left\{ {\boldsymbol{Q}}^0_{ij} \,\, \nabla_i^0\Omega \,\, \nabla_j^0\Omega +\sum_{\alpha=1}^{n} {\boldsymbol{Q}}_{ij}\,\, \nabla_i^{\alpha}\Omega \,\, \nabla_j^{\alpha}\Omega \right\}. \end{equation} \tag{ 66 }$$

These couplings can be incorporated into the isotropic term in equation (58) to yield

$$\begin{equation} \fl H_{\mathrm HR} =\int \mathrm{d}{\hat{\boldsymbol{x}}} \left\{ \frac{1}{2} \boldsymbol{l}^0_{ij}\, \nabla^{0}_i \Omega\, \nabla^{0}_j \Omega + \frac{1}{2} \sum_{\alpha = 1}^n \boldsymbol{l}_{ij}\, \nabla^{\alpha}_i \Omega\, \nabla^{\alpha}_j \Omega + \frac{1}{2} r \, \Omega^2 - \frac{1}{3} w \, \Omega^3 + \cdots \right\}, \end{equation} \tag{ 67 }$$

where

$$\begin{equation} \boldsymbol{l}^0 \equiv \boldsymbol{I} + {\boldsymbol{Q}}^0\quad\mathrm{and} \quad \boldsymbol{l} \equiv \boldsymbol{I} + {\boldsymbol{Q}} \end{equation} \tag{ 68 }$$

are the polymer step-length tensors in the preparation and measurement states, respectively.

The general saddle-point in the presence of a macroscopic deformation gradient Λ, which we denote by $\overline {\Omega }_{\boldsymbol {\Lambda }}({\hat {\boldsymbol {x}}})$, can also be obtained, by minimizing the free energy (67):

$$\begin{equation} \overline{\Omega}_{\boldsymbol{\Lambda}}({\hat{\boldsymbol{x}}}) = q\, \int_z \int_{\tau} \exp\left\{- \frac{\tau}{2} {\boldsymbol{y}}^0 \cdot \boldsymbol{l}^0 \cdot {\boldsymbol{y}}^0 - \frac{\tau}{2}\sum_{\alpha = 1}^n {\boldsymbol{y}}_{\boldsymbol{\Lambda}}^{\alpha} \cdot \boldsymbol{l} \cdot {\boldsymbol{y}}_{\boldsymbol{\Lambda}}^{\alpha}\right\} - \frac{\,\,q}{\quad V^{n}}, \end{equation} \tag{ 69 }$$

where

$$\begin{equation} {\boldsymbol{y}}^0 \equiv {\boldsymbol{x}}^0 - {\boldsymbol{z}} \quad{\rm and}\quad {\boldsymbol{y}}^{\alpha}_{\boldsymbol{\Lambda}} \equiv {\boldsymbol{x}}^{\alpha} - \boldsymbol{\Lambda} \cdot {\boldsymbol{z}}. \end{equation} \tag{ 70 }$$

Note that the Gaussian fluctuations of the particle positions have variance matrices l⁰ in the preparation state and l in the measurement state. The effect of nematic order is therefore to render the thermal position fluctuations of each molecule anisotropic.

Substituting this deformed saddle-point (69) back into equation (67), and carrying out the prescription specified in equation (34a), we find the following elastic free energy for nematic elastomers:

$$\begin{equation} H(\boldsymbol{\Lambda}) = \frac{\mu}{2}\,{\rm{Tr}\,}\,\boldsymbol{l}^0 \boldsymbol{\Lambda}^{\rm{T}} \boldsymbol{l}^{-1} \boldsymbol{\Lambda} \end{equation} \tag{ 71 }$$

which is precisely the neo-classical elasticity theory for nematic elastomers, derived by Warner and Terentjev [3, 9] as a generalization of the classical theory of rubber elasticity.