Generalized Deam-Edwards Approach to the Statistical Mechanics of Randomly Crosslinked Systems

We address the statistical mechanics of randomly and permanently crosslinked networks. We develop a theoretical framework (vulcanization theory) which can be used to systematically analyze the correlation between the statistical properties of random networks and their histories of formation. Generalizing the original idea of Deam and Edwards, we consider an instantaneous crosslinking process, where all crosslinkers (modeled as Gaussian springs) are introduced randomly at once in an equilibrium liquid state, referred to as the preparation state. The probability that two functional sites are crosslinked by a spring exponentially decreases with their distance squared. After formally averaging over network connectivity, we obtained an effective theory with all degrees of freedom replicated 1 + n times. Two thermodynamic ensembles, the preparation ensemble and the measurement ensemble, naturally appear in this theory. The former describes the thermodynamic fluctuations in the state of preparation, while the latter describes the thermodynamic fluctuations in the state of measurement. We classify various correlation functions and discuss their physical significances. In particular, the memory correlation functions characterize how the properties of networks depend on their history of formation, and are the hallmark properties of all randomly crosslinked materials. We clarify the essential difference between our approach and that of Deam-Edwards, discuss the saddle-point order parameters and its physical significance. Finally we also discuss the connection between saddle-point approximation of vulcanization theory, and the classical theory of rubber elasticity as well as the neo-classical theory of nematic elastomers.


I. INTRODUCTION
Randomly crosslinked systems-including rubbers, gels, liquid crystalline elastomers, cytoskeleton networks, and many other related systems-arise widely in nature and are also extensively manufactured by scientists and technologists . Their physical properties are dominated by heterogeneous, system-spanning random networks. A proper specification of this random network structure therefore constitutes the first step towards obtaining a more complete understanding of these complex materials.
In a typical protocol for making randomly crosslinked materials, one starts from a liquid system under particular physical conditions (e.g., of temperature, pressure, solvents, etc.), and randomly introduces crosslinkers that link together nearby particles or polymers. The macroscopic liquid state prior to crosslinking is defined as the preparation state. A realistic crosslinking process necessarily takes some finite interval of time, and is therefore inevitably accompanied by some irreversible responses of the system. Inclusion of these irreversible aspects would make the theoretical modeling of the random crosslinking process extremely complicated. To simplify the analysis and nevertheless capture the essential physics of random crosslinking, we shall therefore follow the original idea of Deam and Edwards [1] by considering an idealized, instantaneous crosslinking scheme in which all crosslinkers are introduced into the system at one particular instant in time. If a sufficient number of crosslinkers is introduced, a system-spanning random network emerges, and endows the material with a non-vanishing shear modulus. A continuous transition such as this, from a liquid to a random solid, is usually called gelation-or vulcanization in settings of pre-existing polymers. Gelation and vulcanization are quite different from usual thermodynamic phase transitions, because they are irreversible and the resulting systems acquire intrinsic spatial heterogeneity. The statistical theory for these transitions is therefore necessarily more complicated.
It is well known that the physical properties of randomly crosslinked materials depend both on the state of measurement and on the state of preparation (i.e, the history of formation). This is, of course, a feature common to all heterogeneous materials, and has been widely appreciated by experimentalists. It was also explicitly pointed out by Deam and Edwards [1] and by de Gennes [2] as early as the 1970's. It therefore follows that the statistical physics of randomly crosslinked materials involves two statistical ensembles: the preparation ensemble and the measurement ensemble. This necessitates a substantial generalization of the standard Gibbs ensemble theory for equilibrium statistical mechanics.
To substantiate this seemingly extreme statement, let us discuss the remarkably distinct properties exhibited by nematic elastomers that have been crosslinked under distinct conditions. Nematic elastomers [3] are crosslinked nematic polymer melts, and their liquid crystalline ordering is strongly coupled to their network elasticity. Two crosslinking protocols have been studied recently by Urayama et al . In one, the system is first crosslinked in the hightemperature, isotropical phase of liquid crystallinity, and is then brought to a lower-temperature, nematic phase. The resulting system (called an isotropic-genesis polydomain (PD) nematic elastomer, or I-PNE) exhibits a polydomain nematic state, with a typical nematic domain size of a few microns. Upon stretching, the nematic domains gradually rotate and deform, and eventually, beyond certain strain value, they align along the external stress so as to form a monodomain (MD) state. The stress response during the PD-MD transition window is extremely soft, being orders of magnitude lower than that of typical elastomeric materials.
In the other protocol, the system is first quenched into the nematic phase, so that the nematic order is well developed locally, but is still frustrated by nematic defects at larger scales (typically beyond 100 µm). Crosslinking carried out in this state freezes in the patterns of nematic defects. The resulting system (known as a nematic-genesis polydomain nematic elastomer, or N-PNE) also exhibits a stress-driven PD-MD transition. The observed stress is, however, within the range of typical elastomer materials, and much higher than that of I-PNEs. The striking differences between I-PNEs and N-PNEs are entirely due to their distinct preparation states. Evidently, a theoretical understanding of these differences demands a statistical theory that involves both the preparation ensemble and the measurement ensemble.
A generic protocol for preparing and measuring a randomly crosslinked system has three steps, each statistical in nature. First, there is a preparation state. We shall denote averaging over the thermal fluctuations in this state by square brackets [ A ] p , where the subscript p stands for "preparation. " Second, there is an instantaneous crosslinking stage, at which crosslinkers are randomly introduced into the system. Averages over realizations of crosslinkers are denoted by an overbar: A . Finally, there is a measurement state. Averages over fluctuations in this state are denoted by angle brackets A m , where the subscript m stands for "measurement. " Note that these three types of average have a causal ordering in time: Any quantity averaged in a measurement state must also be averaged over crosslinker realizations and over the preparation state, as the latter two ensembles constitute the history of formation. By contrast, quantities averaged in the preparation ensemble are not to be averaged over crosslinking, because the network does not exist in the preparation state.
From the above example of nematic elastomers, it can be inferred that the structure of a polymer network in a randomly crosslinked system depends not only on the randomness inherent in the crosslinking process but also on the thermal fluctuations in the state of preparation. Two macroscopically identical crosslinking processes using the identical protocol yield two distinct networks which are only statistically similar. More importantly, two such crosslinking processes carried out under different preparation states yield statistically different network structures, and lead to different properties of the systems, even under the same conditions in the measurement state. It therefore takes both the crosslinking methods and the preparation state to have a full statistical description of the structure of the random network in these materials.
In view of the three steps just mentioned, as well as their statistical nature, there are four types of physical quantities that can be measured in a randomly crosslinked system: (1) Thermodynamic averages in the preparation state: [A] p , [A B] p etc., where the subscripts p stands for "preparation. " (2) Thermodynamic averages in the measurement state: A m p , AB m p etc. Note that these are to be averaged over both their crosslinking randomness and the thermal fluctuations in the preparation state. The purpose of this Paper is to discuss various basic aspects of vulcanization theory, which provides a general methodology for the systematic calculation of all four of these types of diagnostic quantity. The remainder of the Paper is organized as follows. In Sec. II we introduce a toy model of randomly linked anisotropic particles, as a starting point for vulcanization theory. In Sec. III we discuss a generalized Deam-Edwards distribution of network structure, and average the free energy of the network over this distribution. At this stage, we shall see a key theme, viz., that both the preparation and the measurement ensemble appear naturally, coupled to one another, in the effective replicated Hamiltonian. In Sec. IV we discuss in detail various physical quantities that can be calculated using the replicated partition function. We also discuss the significance of the zeroth replica and the connection between replica limit and causality. We also address the distinction between the original Deam-Edwards distribution and our generalization of it. In Sec. V we discuss the physical meaning of the vulcanization order parameter. In Sec. VI we study the mean-field approximation, and from it we derive the classical theory of rubber elasticity, as well as the neoclassical theory elasticity, appropriate for nematic elastomers. Finally in Sec. VII, we draw the concluding remarks by briefly reviewing three basic results about the heterogeneous nature of randomly crosslinked materials, and by indicating some possible future directions.

II. MODEL OF LINKED PARTICLES
To make our discussion concrete, we consider a liquid of particles, each of which may or may not be anisotropic. The particles may, e.g., be entire polymers that we regard in a coarse-grained sense. Under appropriate conditions, the particles may participate in some kind of collective ordering, such as liquid crystalline alignment. Each particle has certain number of functional sites at which the crosslinkers can act. Before crosslinking, the spatial locations of all the functional sites of all the particles, which we index via i, are denoted by c 0 i . We assume that the positions of these sites c 0 ≡ {c 0 i , i = 1, . . . , N } completely specify the microscopic state of the liquid (i.e., there are no relevant degrees of freedom beyond these). The liquid Hamiltonian (normalized by temperature) is denoted by H 0 liq (c 0 ). The corresponding partition function for the liquid system (in the preparation state) is then given (with the shorthand notation   (1,9), (2,8), (3,4), (4,5), (5,9), (6, 7)} The crosslinkers are modeled as Gaussian springs, each connecting a pair of functional sites (i, j) (for i ≤ j to avoid double counting). The network structure is uniquely specified by the list of N (N + 1)/2 "link numbers" k (i,j) that specify the number of springs linking the pair of sites (i, j): χ = {k (i,j) } [14]. Then, the summation over network structures corresponds to summations over all (non-negative) integer values for k (i,j) : Another, equivalent, specification of the network structure is given as follows. First, we specify the total number of springs M . Second, we treat these M springs as if they are distinguishable, and number them using an integer index e = 1, . . . , M . The e th spring can be denoted by an integer pair (i e , j e ) (with i e ≤ j e ) which specifies the pair of functional sites it links together. The network structure is then fully but non-uniquely characterized by an integer M and a list of M integer pairs χ = {(i e , j e )} M e=1 . Because all springs are identical, any arbitrary permutation of the M pairs within the list χ does not change the network structure. We easily find the number of distinct χ that correspond to the same network structure (which is uniquely determined by χ), i.e., the "degeneracy factor": ( The summation over all network structures can then also be represented as a summation over the crosslinker number M , combined with a sum over the list of 2M integers {i e , j e }, Note, however, that the probability of a given network structure P χ must be divided by the degeneracy factor (3) in order to correct for the over-counting. We shall discuss this issue in detail below. For an illustration of a particular network structure, see Fig. 1. For the two equivalent specification of this network structure, see Table II. Let us assert that after crosslinking the system is brought to the measurement state, in which the microscopic locations of the sites are now labeled by c i (rather than cv 0 i of the preparation ensemble). The interaction between pairs of sites linked by springs is modeled via quadratic terms in the separations between sites: This potential energy increases without bound as any two linked sites are separated, and hence it keeps them close to each other at all times. For later purposes, we also define the liquid-state partition function in the measurement state: The partition function and free energy of the crosslinked system are then given, respectively, by both of which depend on the network structure χ. Throughout this paper, we shall set the Boltzmann constant to be unity. Note that, generically, H liq (c) (i.e., the liquid Hamiltonian in the measurement state) differs from its preparation-state counterpart H 0 liq (c). We then average the free energy over realizations of the network structure χ (with probability P χ ) to obtain Assuming that the selfaveraging property holds, this formula also gives the free energy of a typical sample. The calculation (or modeling) of the probability P χ for the network structure χ is the sine qua non for obtaining a proper statistical treatment of randomly crosslinked systems. This we now discuss in detail.

III. GENERALIZED DEAM-EDWARDS DISTRIBUTION OF NETWORK STRUCTURE
A. Statistics of connectivity between a pair of sites 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. Eq. (4)]: In addition, we assume that each functional site c 0 i has an "effective region" of radius , centered around c 0 i . If an end of a spring is introduced at some position x that lies within the "effective region" around c 0 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.
Next, let p k (c 0 i , c 0 j ) be the probability that a given pair of sites {c 0 i , c 0 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: where we have introduced the dimensionless control parameterμ for the crosslinking density [15], and the constant C is determined via overall normalization, see below. Note that p k (c 0 i , c 0 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 2 : The function δ b (x) tend to the Dirac delta function as b → 0. The linking probability (9) can then also be expressed in the following concise form: where we have exchanged the crosslink density control parameterμ for the parameterμ, defined viã It is easy to check that Eq. (11) is properly normalized:

B. Statistics of the connectivity of the network
We assume that the liquid, prior to crosslinking, is in thermal equilibrium with respect to Hamiltonian H 0 liq (c 0 ), where c 0 ≡ {c 0 i } denotes the locations of all N functional sites in the preparation state [16]. The joint probability density function for all N sites (prior to crosslinking) to be at the locations {c 0 1 , . . . , c 0 N } is then given by where the normalizing partition function is given by Now recall that the network structure is characterized by a set of N (N +1)/2 stochastic variables [17], 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 0 prior to crosslinking, these link numbers are fully uncorrelated, and governed by the following conditional probability distribution: where we have used the linking probability given by Eq. (11). Making the further definition we can write this conditional probability as where we have used Eq. (10), Eq. (12), and Eq. (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 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 Furthermore, by making the definitioñ we obtain Hence we see that, apart from a dimensionless factor, P χ is the ratio of two partition functions. Even though the variables c 0 i (associated with preparation state) appear only as dummy variables, the integrations over them in Eq. (20) cannot be carried out explicitly. Below, when we calculate the disorder average of free energy, we shall find that these variables c 0 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.

C. Normalization of probabilities
The conditional probability that there are precisely M linking springs in the system can be obtained by summing Eq. (18) over all possible values taken by the indices (i e , j e ): It is elementary to check that the conditional probability P (M |c 0 ) is properly normalized for any fixed c 0 : To check the normalization of the probability P χ , we sum Eq. (19) over all possible network structures. By Eq. (2), this amounts to summing over all possible linking numbers {k (i,j) }. Further using Eqs. (4) we arrive at The summation and product can readily be calculated, leading to the factor exp −H X (c 0 ), which precisely cancels the corresponding exponential factor. Hence, as expected, we have the normalization condition

D. Averaging over network structures
We are now ready to compute the disorder-averaged free energy, Eq. (7). For this purpose, we use the following identity, which is commonly called the replica trick [18]: We use this to rewrite the average free energy as We calculate the right hand side of Eq. (27) for positive integers n, and subsequently analytically continue n to real values. More specifically, we calculate Z n χ by replicating the variables {c i } a number n times, to obtain an n-tuple (c 1 i , c 2 i , . . . , c n i ). Together with the corresponding variable c 0 i of the preparation ensemble, they form a (1 + n)-tuplê . The shorthandĉ has been extensively used for this (1 + n)-tuple in the vulcanization theory literature. Using Eqs. (6,20), the product of partition functions Z 0,χ (Z χ ) n in Eq. (27) can be written in terms of functional integrals (with the shorthand notation Dĉ ≡ n α=0 Dc α ) as follows: Next, by using the definition Eq. (4) we easily see that where we have used the notation |ĉ α | 2 ≡ n α=0 |c α i | 2 .
We can now substitute Eq. (28) back into Eq. (27), exchange the order in which χ and Dĉ are performed in Eq. (27), and sum over χ for a fixed configuration ofĉ. The summation over the realizations of network structure can be calculated using the following result: where in the last equality, we have used the definitions As we shall eventually take the replica limit, n → 0, the difference between µ andμ [see Eq. (12)] can be ignored.
Combining all these results, we find that the disorder-averaged free energy is given as follows: The Hamiltonian H norm [ĉ] represents an effective short-range attraction between all pairs of functional sites in the replicated coordinate spaceĉ = {c α } n n=0 , resulting from cross-linking. This short-range interaction reduces to a delta function in the limit b → 0.

IV. PHYSICAL QUANTITIES AND CORRELATION FUNCTIONS
We shall use the notation Ψ(ĉ) 1+n to indicate averages of physical quantities with respect to the weight in the replicated partition function Eq. (34b): When there is no risk of confusion, we shall neglect the subscript 1 + n on the average and simply write Ψ(ĉ) . It is understood that at the end of the calculation of such quantities the replica limit n → 0 should always be taken. In the replicated Hamiltonian H 1+n , Eq. (34c), only the last term, H norm (ĉ), couples together various replicas. It can be understood as a short-range attraction between all pairs of the 1 + n replicated positionsĉ = (c 0 , c 1 , . . . , c n ). The strength of this attraction is controlled by the density of crosslinking µ. If µ = 0, i.e., the system is not crosslinked, H norm = 0, and H 1+n reduces to a sum of liquid Hamiltonians, one for each replica: Variables in distinct replicas then become completely decoupled from one another, and the replicated partition function Z 1+n describes the equilibrium liquid in the preparation state together with n copies of equilibrium liquid in the measurement state, with no interactions amongst these liquids.
For systems with crosslinks, µ = 0 and so there are nonvanishing correlations between the variables c 0 and the other c α . These correlations characterize the dependence of the properties of these systems on the history of formation. Let us look systematically at the various physical quantities that can be constructed within the framework of an (1 + n)replica theory, and relate them to the suite of correlators involving the preparation and measurement ensembles discussed in the Introduction. Below, we use the shorthand notation A α ≡ A(c α ) for physical quantities, and we frequently omit the subscript 1 + n on various expectation values, all to ease the notation. 2. Two-replica quantities: A α B β = A m B m (for α = β; α, β = 0) are nonvanishing in the presence of the quenched disorder arising from random crosslinking, and thus give the glassy correlations in the measurement state. The order parameter for the random solid state is an especially important example of such a multi-replica quantity.
Memory correlation functions involving 0-th replica only once, A 0 B α c = A B m p (for α = 0) measure the cross correlation between the preparation state and the measurement state. These memory effects are a hallmark property of randomly crosslinked materials. Vulcanization theory provide the unique framework via which they can be systematically calculated.
In the setting of liquid crystalline elastomers, the memory correlators are the quantities that can best distinguish systems prepared under distinct physical conditions, e.g., isotropic genesis nematic elastomers (IGNEs) and nematic genesis nematic elastomers (NGNEs). Various correlation functions, expressed in both replica notation and physical notation, are summarized in Table IV.

A. Significance of the zeroth replica
We can use Z 1+n with appropriate source terms to calculate the average of a physical quantity A(c 0 ) (≡ A 0 ) that only depends on the degrees of freedom c 0 of the preparation ensemble. As replicas 1, . . . , n play no role in such a II: Various correlation functions in vulcanization theory, showing the correspondence between physical notation and replica notation. (Note that in the table we have 0 = α = β = 0.) In the replica notation, all averages are defined by Eq. (35).

Measurement correlation
Glassy correlation Memory correlation 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 . We also have to delete the integrals over c α (for α = 1, . . . , n) in the functional integrals (34b) and (35). Therefore, we have As a result, A 0 is indeed the average in the preparation state, as expected.

B. Replica limit and causality
In the replicated partition function, Eq. (34b), there is one copy of the preparation ensemble (the 0-th replica), and n copies of the measurement ensemble. In Fig. 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 Eq. (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 0 in fluid Hamiltonian H 0 liq of the preparation state [see Eq. (34c)], there will be changes in the physical quantities in the measurement state: This is, of course, entirely reasonable, as the properties of polymer networks do depend on their history of formation. 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.
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. 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 Substituting this back into the expectation value Eq. (36), taking the derivative with respect to J, and finally setting J = 0, we obtain where in the last equality we have assumed that permutation symmetry amongst replicas 1 to n remains intact. The right hand side of Eq. (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 Eq. (39) would even be proportional to n, and hence would vanish in the replica limit.

C. Deam-Edwards revisited
In the replicated Hamiltonian (34c), the parameters of the preparation state only appear in H 0 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 0 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 0 liq (c 0 ) − H norm (c 0 ) instead of just H 0 liq (c 0 ): Consequently, the replica Hamiltonian H 1+n is not invariant under permutations of the replicas that mix the 0-th 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: This is completely reasonable, as the random network structure does affect physical quantities in the measurement state but cannot not 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 The corresponding (1 + n)-replica Hamiltonian is then given by 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.

D. Comparison with spin glasses and other disordered systems
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]: 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: 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 history of formation. Hence, the proper statistical modeling of such systems requires a separate statistical ensemble.

V. ORDER PARAMETER FOR VULCANIZATION
The crosslinking-induced interaction (34d) couples distinct functional sites. It can be decoupled using the standard Hubbard-Stratonovich transformation. To proceed, we first notice that Gaussians are closed under convolution. Therefore, the following identity about the (1 + n)d-dimensional soft delta function δ b (x) [defined in Eqs. (10,33)] can readily be established: Equation (34d) can then be written as is the microscopic collective density of functional sites, coarse grained by the Gaussian wave packet δ b (x). It depends explicitly on the microscopic configuration of the replicated system. Our definition of the collective field differs from those in the other literature on vulcanization theory [6] via a multiplicative factor of N . We can now introduce an auxiliary field Ω(x) ≡ Ω(x 0 , . . . , x n ) via the following Hubbard-Stratonovich transformation: where C is a trivial constant that will be neglected henceforth. Insertion of this result into Eq. (34b) leads to Let us further define a modified liquid state partition function viã and the average · · · 0 1+n as Then we can rewrite Eq. (45) as where the effective Hamiltonian H eff (Ω) is defined as We can now define the average of functions of Ω using the equilibrium partition function associated with H eff (Ω): The following identity is a well-known property of the Hubbard-Stratonovich transformation: which shows the physical significance of the collective field Ω(x), namely that its average is the joint probability density function of the following 1 + n events: in the preparation state there is a particle at some point x 0 , whereas in n independent measurements on the measurement state the same particle is found at positions x 1 , . . . , x n , respectively. Using the physical notation discussed in the Introduction, this probability density can be expressed as where, as usual, the subscripts p and m respectively denote the preparation and measurement states, and the overline indicates an average over realizations of crosslinkers. Ω(x) encodes information about the localization of particles in the gel phase, and consequently is called the vulcanization order parameter. The definition of Ω(x) in this work differs from that in some of the other vulcanization literature via a trivial additive constant and a multiplicative constant N .

A. One replica sector and marginal distributions
We can integrate out the n variables x 1 , . . . , x n in the order parameter (51), leaving only one variable x 0 , and thereby obtain a one-replica distribution 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 δΩ 0 (x 0 ) = Ω 0 (x 0 ) − Ω 0 (x 0 ) 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: 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 Eq. (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 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 Ω(x).

VI. LANDAU THEORY AND SADDLE-POINT APPROXIMATION
Eqs. (48, 50) provide a formal approach to calculate all relevant physical quantities. Our problem is therefore reduced to the calculation of a functional integral over Ω(x). At the saddle-point level we approximate the disorderaveraged free energy (34a) using the minimum of the effective Hamiltonian H eff [Ω], Eq. (49). Even after making this approximation, the problem remains difficult, because we do not know the effective Hamiltonian (49) exactly, let alone its minimum. We therefore proceed to expand the exponential in Eq. (49) in a Taylor series in Ω. At this stage, it is convenient to shift the order parameter (51) by a constant N/V 1+n (where N is the total number of functional sites), so that the it satisfies the condition dx Ω(x) = 0.
Both in the liquid phase and in the gel phase, the system is uniform in density. Therefore the saddle-point order parameter satisfie: We therefore only have to consider the subspace of Ω(x) in which Eq. (57) is obeyed. This subspace is usually called the higher replica sector.
The form of Landau free energy given in Eq. (49) is generic for all isotropic networks: The coefficients in the Taylor expansion do of course depend on the short-scale structure of the constituent particles. It may appear surprising that the cubic term in Eq. (58) has a negative coefficient, which seems to suggest that the Landau theory is unstable, at lest for uniform values of the order parameter. By Eqs. (53,54), however, a large and uniform Ω(x) also implies the large deviation of the particle densities from their equilibrium values in each replica [and therefore would violate Eq. (57)], and this would incur a large free-energy penalty. The Landau theory of vulcanization is therefore locally stable. Minimization of Eq. (58) leads to the following saddle-point equation: This equation should be solved subject to the constraint (57), and is therefore necessarily nonuniform in replicated space.
A. Order parameter at the saddle-point 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 Ω(x). For r > 0, the saddle-point value of Ω is trivially zero, corresponding to the liquid state. For r < 0, the trivial saddle-point becomes unstable, and a nontrivial (i.e., replica-space dependent) saddle-point value emerges, corresponding to the amorphous solid state. It has the following form of superposed Gaussians [19]: 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 Refs. [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 Eq. (60): The vector z can be interpreted as the mean position of a particle in the gel fraction, and is localized around x 0 , 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 Eq. (51) is nonzero. The order parameter Ω(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 Eq. (52) combine to give 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 0 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 0 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 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 0 is uniformly distributed over space, and therefore for an arbitrary function We now substitute Eq. (62) back into Eq. (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 Eq. (60).

B. Translational symmetry and order parameter
The 1 + n replicated Hamiltonian H 1+n , Eq. (34c), possesses the symmetry of independent translations of the replicas, i.e., it is invariant under the α-dependent translations: This symmetry is also possessed by the average order parameter Ω(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, Eq. (60). Translations of all particles in one replica leads to a distinct but energetically degenerate order parameter: Ω(x 0 , . . . , x α + u, . . . x n ) = Ω(x 0 , . . . , x α , . . . x n ).
This reduction in 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 Eq. (60): Ω(x 0 , x 1 , . . . , x n ) = Ω(x 0 + u, x 1 + u, . . . , c n + u).
This macroscopic translational invariance (MTI) [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.

C. Classical theory of the elasticity of rubber
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 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 VII. CONCLUDING REMARKS: WHY VULCANIZATION THEORY?
It is quite satisfactory that the classical theory of rubber elasticity and the neo-classical elasticity theory of nematic elastomers can be derived from the vulcanization theory as saddle-point approximations. These theories were originally established by studying the statistics of single polymers inside a network, with two key assumptions [3,8]: (1) that strain deformations are affine at all scales; and (2) that the polymer statistics are Gaussian. These theories have attained classical status because of their remarkable successes in explaining the basic properties of isotropic and nematic elastomers. Nevertheless, the starting point of single polymer statistics substantially restricts the generalizability of these theories. For example, it is hard to see how one might systematically explore the nature of heterogeneities in random polymer networks using these molecular-level approaches.
Vulcanization theory was developed to understand the heterogeneities of randomly crosslinked systems. In the face of its (perhaps heavy!) machinery-such as field theory, the replica technique, and multiple statistical ensembles-if all that could be obtained from it were re-derivations of classical results then it would be hard to argue that vulcanization theory has been worth investing in. Fortunately, however, several fundamental new results have already been obtained using vulcanization theory, concerning the heterogeneity of rubbery materials: (1) the distribution of localization lengths p(τ ) in isotropic networks [10]; (2) the distribution of internal random stresses in isotropic networks [11,12]; and (3) the statistics of random fields and memory effects in isotropic-genesis nematic elastomers [13]. We refer the reader to the literature for further details about these topics. We think it is fair to say that we are still in an early stage in the unraveling of the statistical physics of network heterogeneities and memory effects in randomly crosslinked materials. We hope that fresh and interesting results will be available to report in the near future.
We thank our numerous co-workers and several constructive anonymous referees for valuable collaborations and guidance. This work was supported by NSFC (China) via grant nos. 11174196 and 91130012, by the U.S. NSF via DMR 09-06780 and DMR 12 07026, and by the Institute for Complex Adaptive Matter.