Theory of Flexible Polymer Networks: Elasticity and Heterogeneities

A review of the main elasticity models of flexible polymer networks is presented. Classical models of phantom networks suggest that the networks have a tree-like structure. The conformations of their strands are described by the model of a combined chain, which consists of the network strand and two virtual chains attached to its ends. The distribution of lengths of virtual chains in real polydisperse networks is calculated using the results of the presented replica model of polymer networks. This model describes actual networks having strongly overlapping and interconnected loops of finite sizes. The conformations of their strands are characterized by the generalized combined chain model. The model of a sliding tube is represented, which describes the general anisotropic deformations of an entangled network in the melt. I propose a generalization of this model to describe the crossover between the entangled and phantom regimes of a swollen network. The obtained dependence of the Mooney-Rivlin parameters C1 and C2 on the polymer volume fraction is in agreement with experiments. The main results of the theory of heterogeneities in polymer networks are also discussed.


Introduction
Polymer networks and gels belong to a unique class of materials that have the properties of both solids and liquids. Such "soft solids" are elastically deformed at macroscopic scales, while at short distances and short times the network strands experience liquid fluctuations, which are responsible for the exceptional properties of polymer networks, such as the reversibility of huge deformations (when stretched up to 3000%, see [1]).
In this work, I review the main microscopic (molecular) approaches to the description of the elasticity of flexible polymer networks and establish relationships between these approaches. Phantom networks are considered in Section 2, and networks with topological entanglements are studied in Section 3. Spatial heterogeneities developed in swollen and deformed polymer networks are investigated in Section 4.
The classical theory proposed by Flory more than half a century ago, Ref. [2] suggests that polymer networks have a tree-like structure. Although this theory takes into account the presence of loops in the network, it is assumed that they all are infinite in size. Such an approach is insufficient to describe the thermodynamics of polymer networks with loops of finite size, and explicit consideration of these loops is required. The mathematical model of such networks is based on the replica method. This approach takes into consideration that the properties of the network depend not only on the "current" conditions in which the network is placed but also on the conditions of its preparation since it is these conditions that determine the molecular structure of the formed network. The replica approach takes into account the excluded volume interaction of linear chains both at preparation and experiment conditions by generalizing the theory of the "n → 0-component" field ϕ 4 proposed by de Gennes.
The main properties of the order parameter ϕ introduced for the description of such "soft solids" are discussed in Section 2.
The conformations of the chains in the network are described by the combined chain model [3]. This model is based on the concept of virtual chains, which determine the effective elasticity of the tree-like structures in the network. In polydisperse networks, the distribution of virtual chains is also polydisperse. I show that the order parameter in the replica network model stores information about this distribution. The found solution of the replica model is used to calculate the distribution function of the lengths of virtual chains in polydisperse networks.
The structure of real networks differs significantly from the classical picture of ideal trees. Typical loops of such networks have finite dimensions and strongly overlap with each other. The impact of such loops on the conformations of the network chain is described by the extended combined chain model [4]. I discuss the connection of this model with the predictions of the replica theory, which offers an analytical description of such mutually overlapping loops. Along with typical loops, the network may contain topological defects, which are also taken into account by the replica method. The condition for the loss of elasticity by polymer networks, both due to the presence of topological defects and near the gel point is discussed.
With an increase in the length of polymer chains, the effects of topological entanglements become significant. The properties of entangled polymer liquids are significantly different from the properties of entangled soft solids-polymer networks. While non-concatenated rings in the melt are compacted into fractal loopy globules [5], internal stresses in stretched polymer networks prevent such a collapse. Entangled networks are described by the non-affine tube model that generalizes the concept of virtual chains to the case of entangled polymer systems [6]. I discuss the physics of deformation of network strands in this model in Section 3.
Anisotropic deformations of the entangled network lead to slippage of chains along the contour of the entangled tube, which is described by the slip tube model [3]. This model suggests the additivity of the phantom and entangled contributions to the free energy of the network. Although this additivity approximation describes well the uniaxial deformations of the network in the melt, it cannot be used to describe the crossover of the Mooney Rivlin dependence between entangled and phantom deformation regimes observed during swelling of the polymer network. I propose the generalization of the slip tube model to account for this crossover and show that the obtained concentration dependences of the Mooney-Rivlin parameters C 1 and C 2 are in agreement with the experimental data.
The irregularity of the molecular structure of polymer networks leads to the presence of spatial inhomogeneities. The main results of the theory of heterogeneities in polymer networks are discussed in Section 4.

Phantom Networks
The polymer network is obtained by crosslinking linear chains. In models of phantom polymer networks, it is assumed that in the process of thermodynamic fluctuations, polymer chains can freely "pass" through each other. The phantom approximation works well only for not too long network strands; otherwise, topological restrictions related to the mutual impermeability of polymer chains should be taken into account. In this section, I discuss the main approaches to the theoretical description of phantom networks. The elastic free energy of a network of Gaussian chains is quadratic in the deformation ratios λ α = L α /L 0α of the network dimensions along the principal axes of deformation α = x, y, z in the deformed state (L α ) and under conditions of its preparation (L 0α ): The interaction of strand monomers is usually taken into account by imposing the incompressibility condition λ x λ y λ z = 1/φ. (2) where φ is the polymer volume fraction. In the case of uniaxial deformation we have In the affine network model, it is assumed that the ends of each network strand are pinned to a "non-fluctuating elastic background" which deforms affinely with the network surface [3]. The elastic modulus of such a network is proportional to the density ν of its strands, and kT is the thermal energy. Equation (4) is attractive for its simplicity. However, in real networks, the strand ends are not fixed in space-they are connected to neighboring strands through cross-links. The elastic modulus G depends on the molecular structure of the network and the polydispersity of its strands. At first glance, the "bookkeeping" of the order of connection of all network strands with each other and the numbers of monomers of each of these strands seems impossible for networks of macroscopic sizes. However, this is precisely what the replica method does.

Replica Method
The replica method was proposed a long time ago [7] and was used later to calculate correlation functions of density fluctuations in polymer networks [8]. The main goal of this approach is to build a mathematical model that takes into account most of effects inherent to real polymer networks. Unlike widespread numerical simulations, which can also take these effects into account, the replica method is an analytical theory, which is more amenable for understanding network elasticity as a function of many different molecular features.
When the polymer network is deformed, its structure remains unchanged, formed under the conditions of its preparation. The main idea of the replica approach is to consider the polymer network in an expanded space with coordinatesx = (x (0) , x (1) , . . . , x (m) ) of dimension 3(1 + m). The replica space consists of an "initial system" under conditions of the network preparation (with coordinates x (0) ) and m replicas of the "final system" under experimental conditions (with coordinates x (k) , k = 1, . . . , m), see Figure 1. The replica system consists of the initial system (k = 0)-the polymer network at preparation conditions (blue cube) and m identical replicas of the final system (k = 1, . . . , m)-the network, deformed by factors λ α along the principal axes α = x, y, z of deformation (yellow cuboids). The polymer network in the replica space of dimension 3(1 + m) is mainly localized inside the (green) cylinder with the diameter R aN 1/2 directed along unit vectors {ê α }, corresponding to affine deformation of the polymer. The networks (and all their monomers and strands) in the initial system and any of the m replicas of the final system are the projections of the network in the replica space onto the corresponding subspaces k = 0, . . . , m.
In each of these systems, the polymer network with all of its monomers and bonds is the projection of the network in the replica space onto the corresponding subspaces. Therefore, in all these systems, the network structure is exactly the same, while the conformations of network strands can vary significantly. The free energy F of the final system (a deformed network) is expressed through the analytic continuation to m = 0 of the free energy F m of the replica system [8,9]:

Liquid-Solid Order Parameter
According to de Gennes, polymer networks are "soft solids", which emphasizes the presence of liquid and solid-state degrees of freedom interacting with each other. In the Landau approach, solids are described by an order parameter-a number (in general case, a tensor), equal to zero for a liquid and nonzero only in the solid phase. Unlike ordinary low molecular weight solids, the properties of polymer networks depend on the characteristics of both the initial and final systems. Therefore, the order parameter for polymer networks is not a number, but a function of the difference α of final and initial systems. This dependence reflects the dual nature of polymer networks: At scales larger than the characteristic size of the network cycles, the network deforms affinely as an ordinary solid, see Figure 1. On a smaller spatial scale, the network chains fluctuate in space and have conformations similar to the chains in a polymer liquid. The unique properties of polymer networks are determined by the interaction of their solid-state and liquid degrees of freedom [10].
In the liquid phase, the order parameter is independent of coordinates. In the case of Gaussian polymer networks, the order parameter ϕ (ς) is a function of a single variable ς = 1 2 whereê α is a unit vector along the direction x α of the affine deformation in the replica space.
α ), the order parameter ϕ (0) is determined only by the conditions of the network preparation and does not depend on the experimental conditions. This last dependence is encrypted in the dimensionless order parameter χ (t) = ϕ (ς) /ϕ (0) which is the function of the dimensionless variable t = ς/R. In the case of well-developed networks obtained by crosslinking polymer chains with f -functional monomers far beyond the gel point, the characteristic radius is R bN 1/2 , see Figure 1. Here b is the monomer size andN is the average number of network strand monomers.
In the case of monodisperse networks, the function χ (t) was calculated by Edwards (see Equation (3.19) of Ref. [11]): In the case of polydisperse networks, this function χ (t) is determined by the differential equation For large t 1, the function χ (t) decreases more slowly compared to the monodisperse case (7), as a stretched exponent For small t, the solution of Equation (8) can be approximated by a power function Substituting it in Equation (8) and equating to zero the first two expansion coefficients (t and t 2 ) in powers of t, we find The replica method is also called "replica trick" because of the hidden meaning of its mathematical constructions. As will be shown below, Equation (8) describes the relationship between distribution functions of the molecular trees that characterize the structure of the polymer network. Note that Equation (8) is obtained by estimating the Hamiltonian (A2) in the Appendix A using the steepest descent method. Therefore, it corresponds to the mean-field approximation, which works due to the presence of a large number of overlapping network strands.

Overlap Parameter
Like in liquid polymer systems, in a typical network there are many overlapping network strands. An important characteristic of the network is the overlap parameter, which is defined as the number of network strands within the volume R 3 pervaded by one network strand. In the case of a solution of chains with monomer density ρ (0) , P (0) ρ (0) R 3 /N. In Gaussian networks, the size of the network strand with Kuhn length b is R bN 1/2 , and with a large parameter P (0) , polymer networks obtained far from the gel point can be considered as P (0) overlapping yet independent "elementary" networks, the strands of which consist of N monomers. There is on the order of one network strand per volume R 3 in the "elementary" network, and its elastic modulus is estimated as G elem kT/R 3 , kT is the thermal energy. Since the real network consists of P (0) overlapping elementary networks, its elastic modulus is P (0) -times larger [12]: The polymer network can also be represented as P (0) polymer trees (or "layers") overlapping with each other. The perfect network model assumes an infinite number of layers, and also that these trees have an infinite number of generations. In real networks P (0) < P e is finite and not too large (the entanglement overlap parameter P e ∼ 20-30, see Section 3 below), and the different root trees are interconnected by loops. The number of network strands forming a closed loop of minimum length, which binds together different layers of the network, depends logarithmically on the overlap parameter [4], The case P (0) ≈ 1 with the loop size l ≈ 1 describes the "elementary networks" in Equation (13). In accordance with the results of numerical simulations [13], the loop length l decreases with dilution and with increasing functionality f of cross-links. Note that the minimal size loops do not directly determine the network elasticity since each elastically effective chain is simultaneously part of a large number of loops. The cumulative effect of these typical loops of the network is self-averaged, and can be described in the effective mean-field approximation, see Equation (37) below.
The overlap parameter also determines the density fluctuation of the network obtained in the case of incomplete conversion above the gel point, at p > p c . In the mean-field approximation, the density of gel monomers is and its connectivity radius is estimated as ξ b [N/ (p − p c )] 1/2 . Fluctuations in the number of monomers of the polymer network on the scale of the connectivity radius occur by attaching or detaching sol molecules with a characteristic number of monomers Therefore, relative fluctuations in the gel density are estimated as In the case of a large overlap parameter, these fluctuations are small outside the narrow Ginzburg region,

Combined Chain Model: Monodisperse Networks
The elasticity of the polymer network has an entropic origin and is determined by the change in strand conformations during the network deformation. The conformations of a strand depend on its interaction with the elastic environment, which can be described by virtual chains attached to the ends of this strand. The network strand together with the two virtual chains is called a combined chain. The ends of this combined chain are attached to an affine deformable non-fluctuating background, see Figure 2b In the perfect network model, it is assumed that the network has a tree structure, whereas all its cycles are of infinite size, see Figure 2a. Each of the trees attached to the ends of the network strand is modeled by a virtual chain. The effective number n of its monomers is related to the numbers n i of virtual chain monomers on the next generation of the tree i = 1, . . . , f − 1 and the corresponding numbers of strand monomers N i as In the case of a monodisperse network with a fixed number N i = N of monomers of its strands, all n i = n and the solution of Equation (19) determines the number of monomers of the virtual chain [14] The average vector R between the ends of the network strand is determined from the force balance condition. In the case of a monodisperse network, where n is the number of monomers of the virtual chain, Equation (20). The mean square fluctuation of the vector R is determined by the parallel connection of the real chain and two serially-connected virtual chains The vector X between the ends of the combined chain is deformed affinely with macroscopic network deformation. The elastic modulus of the perfect network with monomer density ρ (0) and chain concentration ρ (0) /N is Here ν = ρ (0) /N is the concentration of strands and (1 + 2n 1 /N) −1 is the strand fraction of the combined chain. A comparison of expressions (23) and (4) shows that the elastic modulus of the perfect network is less than the modulus of the affine network by a factor of 1 − 2/ f . It is generally accepted that the 2/ f factor is associated with the presence of fluctuations in the strand size, which are also proportional to this factor (see Equation (22)). In fact, fluctuations have nothing to do with it [15]: an exact elastic modulus of Gaussian phantom network is different from that of the perfect network model and coincides with the elastic modulus of the classical non-fluctuating grid, in which each strand is replaced by a corresponding elastic thread with the same elastic stiffness coefficient [16]. This coincidence is due to affine deformation of the average distances between the cross-links in such networks.
Equation (23) can be recast in more general form with the network modulus proportional to the difference of the number densities of network strands ν and cross-links µ = 2ν/ f , since there are f /2 strands per crosslink [14]: This expression is quite universal; with general densities ν and µ, it describes polydisperse networks, as well as networks with incomplete conversion, prepared above the gel point far from the Ginzburg region defined by Equation (18). In general, µ − ν makes a sense of the cyclic rank of the network of unit volume. To define it, mentally cut one of the strands so that the network does not break into two disconnected parts. The cyclic rank is equal to the maximum number of such cuts, and it makes sense the number of independent network loops. An important limitation of expression (24) is the assumption of a model of perfect networks about the infinite sizes of all of the network cycles. In real networks, not all loops transmit stress in the network. For example, primary loops, connected to the network at only one crosslink, are not capable of permanently supporting a stress. To exclude the contribution of such elastically ineffective loops and dangling chain ends, it was proposed to account in Equation (24) only for the elastically effective strands and crosslinks. Elastically effective strands deform and store elastic energy upon network deformation [17]. Elastically effective crosslinks connect at least two elastically effective strands [18]. The elastic modulus of real networks can significantly differ from such a modified expression (24), since each of finite size loops is characterized by its "elastic effectiveness", which depends on the type of the loop. The replica method provides a universal tool for accounting for such effects, see Equation (37) below.

Combined Chain Model: Polydisperse Networks
Real networks are polydisperse, and therefore, the numbers of monomers n of virtual chains are random variables. Information on the distribution of the lengths of virtual chains is important for studying the distribution of tension in real chains of the network. Below I calculate the distribution function p (n) of virtual chain lengths. The distribution of inverse variables s = 1/n is described by the function Appendix B provides an algorithm for calculating this distribution for an arbitrary distribution P (N) of the number of monomers N of the real network strands. The problem is reduced to solving a system of nonlinear integral equations. In the most interesting case of the exponential distribution P (N) = e −N/N /N, the asymptotic behavior n N of the solution of these equations can be found Of course, knowledge of only asymptotic is not enough to describe conformations of a chain in the polydisperse network. The replica method offers a much more constructive approach, allowing us to calculate this function over the entire range of n values. One can show, generalyzing calculations of Ref. [9] that the distribution function q (s) in Equation (25) is determined by inverse Laplace transform of the function χ f −1 (t) with the dimensionless order parameter Due to the boundary condition χ (0) = 1 for the function χ (t), the distribution function q (s) is normalized to unity. To better understand the meaning of Equation (27), we substitute the exponential function χ (t) for the monodisperse networks, Equation (7), into this equation. The result is a δ-functional distribution, in accordance with expression (20) for the number of monomers n = 1/s of virtual chains.
In the case of polydisperse networks, this correspondence, Equation (27), establishes a connection between the standard method of distribution functions and the replica approach. Since χ f −1 (t) is the Laplace transform of the convolution functionχ ( f −1) (t) (see Appendix B), the function χ (t) is the Laplace transform of the distribution functionχ (s ) defined by Equation (A7). Equation (8) for this function corresponds to Equation (A8) of the replica approach for the exponential distribution of the numbers of monomers of real strands.
Substituting the asymptotic expression (9) for the function χ (t), into Equation (27), we find a more accurate asymptotic expression for the distribution function Using the function χ (t) from Equation (10), we get Thus, the distribution p (n) of virtual chains vanishes with all its derivatives as n → 0 and decreases exponentially on the scaleN/ ( f − 2) 2 for large n. Calculating the average number of monomers of the virtual chain with the distribution (30), we find Note that the obtained distribution of the lengths of the virtual chains is narrower than the initial distribution of the lengths of the network strands.

Generalized Combined Chain Model
A perfect network assumes a tree-like structure on all spatial scales. In real networks, due to the excluded volume effect, trees cannot "grow" to infinity (the Malthus effect), since too large trees cannot fit in real 3D space. Therefore, in real networks, the size of a typical tree is finite, and the network consists of a large number of loops of finite length, see Figure 2c. They are strongly overlapped and interconnected with each other. The impact of such typical loops, which are responsible for the solid-state elasticity of the network, can be described by the generalized combined chain model.
In this model, in addition to two virtual chains of n 1 monomers, the combined chain includes an additional virtual chain of n 2 monomers, see Figure 2d. This virtual chain represents an effective elastic of all loops of finite size in the network, shunting the real strand of N monomers. The average end-to-end vector R of the network strand is related to the end-to-end vector X of the combined chain via the force balance condition, generalizing Equation (21): In this model, the mean square fluctuation of the vector R is The second equality can be treated as the contribution of two effective virtual chains connecting the network strand to the non-fluctuating background. They have renormalized number of monomers corresponding to parallel connection of two effective chains: with 2n 1 and n 2 monomers. The elastic modulus of this model is The factor (1 + N/n 2 ) −1 is the fraction of energy related to the network strand when it is connected in parallel with the effective chain of n 2 monomers, see Figure 2d. This result can be rewritten in the form where 2n is the number of monomers of the effective virtual chains, Equation (34). The first term in this expression has the same meaning as for the perfect network in Equation (23). The negative second contribution in Equation (36) is due to the finite size of typical loops in the network. The number of monomers of the shunt virtual chain, n 2 , can be calculated using the replica method.

Finite-Size Loops of Real Networks
In real networks, there are always cyclic defects of finite dimensions. Sparse structural defects can be taken into account within the framework of the perfect network model, assuming that they are connected with the affinely deformed non-fluctuating background through the root trees. These trees can be modeled by the virtual chains [3]. Since part of the network strands is "spent" on creating the cyclic fragment itself, the length of the effective chain linking this fragment to the elastic non-fluctuating background is different from Equation (20). In polydisperse networks, the distribution function of such virtual chains is determined by Equation (27), in which f − 1 has the meaning of the number of branches on the first generation of this tree. In the ideal defect gas approximation, the contribution of the structural defects to the elastic modulus of a perfect network was calculated in works [19][20][21][22]. Note that typical network loops are ignored in this approximation, which takes into account only explicitly treated loops of small concentration.
In real networks, there are both loops, binding different network layers (see Equation (14)) and topological defects, which are not sparsely distributed. The replica method allows calculating the effect of both typical loops and cyclic defects of arbitrary concentration on the elasticity of the network. The elastic modulus of the network obtained by random end-linking polydisperse chains by cross-links with functionality f and concentration ρ The functionality f i of the cyclic defect defined as the number of "external" network strands joined to the cyclic fragment. For the ring fragment with i "internal" strands and cross-links, , since two of the functional groups of each cross-link are involved in creating the ring. In the case of incomplete conversion, p < 1 (but still far from the gelation threshold, |1 − p| 1), in expression (37), instead of functionalities f and f i , their average values should be substituted, p f and p f i . Such a dependence of the network elastic modulus on the conversion is in agreement with numerical simulations of Gaussian networks, performed in Ref. [16].
The fractions x i of cyclic defects (small fragments of the network with i 1 "internal" cross-links) depend on the type of reactions leading to the formation of the network [23,24]. All x i universally depend on one dimensionless parameter x 1 characterizing the conditions for network preparation [25]. The primary loops with concentration x 1 are tied to the network at only one crosslink. The denominator of the first term in expression (37) describes an increase in the effective number of monomers between cross-links due to primary loops,N →N (1 + x 1 ). Therefore, although primary loops cannot bear shear stress in the network, they renormalize the length of the elastically effective chains, which deform and store elastic energy upon network deformation [7,9].
For general f > 2, the first term in Equation (37) can be interpreted as the contribution of elastically effective network strands with renormalized number of monomers due to the presence of topological defects in the network-primary loops and cyclic fragments of finite size. Only this contribution is predicted in the classical model of phantom networks. Its expansion in a series in the parameters x i with accuracy up to first order terms reproduces the result of the so-called "network theory with strand pre-strain" [21]. At network preparation conditions, the strands of a cyclic fragment are contracted and the surrounding strands are stretched (pre-strained) with respect to Gaussian sizes of these chains [22]. These effects are "automatically" taken into account in the framework of the replica approach in Equation (37).
The last term in Equation (37) is always negative and gives an impact of a large number of typical loops of the network. The factor Q (0) in the second term of expression (37) determines the probability of the formation of typical closed loops. These loops shunt elastically effective chains of the network, an effect that is taken into account by an additional effective chain in the model of the generalized combined chain in Figure 2d. Unlike perfect networks, the monomers of real networks mutually repel each other, "crowding out" the network loops on small length scales, see Figure 2. Therefore, the elasticity of a polymer network, Equation (37), explicitly depends on excluded volume parameter v (0) , characterizing monomer interaction at the preparation conditions. This effect is ignored by the classical theory of polymer networks.
Note that the elastic modulus of the polymer network G in Equation (37) vanishes at a finite fraction x 1 of cyclic fragments. Such a network with overlap parameter P (0) 1 has tree-like connections of overlapped cyclic fragments with a finite monomer density. To better understand this result, we compare it with the more familiar case of a disordered low-molecular-weight solid. In this case, the network fragments cannot overlap since the overlap parameter P (0) = 1. Therefore, the elastic modulus turns to zero only at the point of percolation transition, at which the network density also vanishes.
A similar effect was observed in numerical simulations of networks obtained near the gel point, which is an analog of the percolation transition in low-molecular-weight solids. It is shown that the elastic modulus of such networks vanishes at a finite density of the network monomers [26]. One can estimate the corresponding shift in critical conversion p cµ at which G = 0, compared with the gel point at the conversion p cγ , at which the average degree of polymerization of soluble molecules diverges. For end-linked networks, the looping probability is Q (0) (bN 1/2 ) −3 ∼ 1/(NP (0) ). Dropping all the numerical factors, we find from Equation (37) the mean-field estimate for the elastic modulus Here L N/ (p − p cγ ) 2 and ρ is the density of gel monomers, see Equations (16) and (15). Therefore, G = 0 at p cµ − p cγ ∼ 1/ √ P (0) , in accordance with the result of numerical simulations [26]. A more rigorous description of this effect requires the inclusion of fluctuation corrections. Inside the strongly fluctuating Ginzburg region, Equation (18), the overlap parameter of elastically effective chains is small (which is why strong density fluctuations develop in this region), and they can be described by critical exponents of percolation theory, see Reference [26].

Entangled Networks
The physics of topological entanglements is perhaps one of the most "entangled" sections of physics of polymer networks since it is impossible to give a sufficiently rigorous and constructive formalism that adequately describes entanglements for networks of macroscopic size. Therefore, the consideration of topological entanglements is usually carried out using some uncontrolled assumptions, allowing to switch from substantially "multi-chain" problem to considering the conformations of one network strand. The transfer of elastic shear stresses from individual strands of a network to its solid-state degrees of freedom (represented by the affinely deformed non-fluctuating background) is described by virtual chains. In phantom networks, such a transfer is performed only through cross-links at the strand ends. In entangled networks, the stress is carried by all entangled segments of the strand due to its topological interaction with the loops formed by adjacent network strands.
Entangled segments are characterized by the large overlap parameter where N e is the number of monomers of the entangled strand. In the case of flexible polymers, the parameter P e 20-30. The number of monomers in an entanglement strand N e and the overlap parameter P e depend on the polymer chemistry [27].
As shown in [5], in the melt of non-concatenated rings the condition P e = const is satisfied at all scales down to the distance a 0 bN 1/2 e between adjacent entanglements. Entangled loops of different sizes overlap with similar size loops at the same overlap parameter P e . As a result, ring molecules are packaged into the self-similar structures called fractal loopy globules. Unlike ordinary globules, in which chains with Gaussian statistics experience random reflections from the surface of the globule, the ring chains in the fractal loopy globules are not Gaussian and form loops at all scales, starting from a size of the entangled segment. Ring sections with the number of monomers g < N e smaller than entanglement strands are still Gaussian with size r bg 1/2 . Larger subsections of a ring are compressed with fractal dimension d f = 3 and size up to the size R of a loopy globule in a melt Large typical values of the parameter N e 1 allow developing the mean-field theory of entanglements in polymer networks, since the total effect of large number P e 1 of loops interacting with the entanglement segment is effectively averaged. The strands of the entangled polymer networks are confined in an effective entangled tube with a diameter a, which, under the conditions of the network preparation, is equal to the size a 0 of the entangled segment. The tube rotates by a random angle of about 90 • on the scale a 0 , see Figure 3a.
Here r α (s) = x α (s) /λ α and X α (s) (α = x, y, z) are undeformed coordinates of two ends of the s-th virtual chain. As shown in [3], with this choice of the potential of virtual chains, they do not contribute to the elastic stress of the polymer network, which is determined only by the contribution of all real polymer chains.

Physics of Entangled Network Deformation
Consider the case of a highly entangled network, N/N e 1, with many entanglements per network chain. As shown in [28], the network deforms affinely like a solid only at spatial scales exceeding the affine length R a f f α , see Figure 3. On shorter scales, chains take on liquid-like conformations. In the case of anisotropic deformation of the gel by the ratios λ α along the axes α = x, y, z, the affine length is also anisotropic and depends on the direction α: We define an affine strand of the size R The solution of Equations (45) and (46) is Thus, the entanglement efficiency decreases with a network stretching, and the entanglements are not effective at N a f f α > N; such an entangled network behaves like a phantom one. The free energy of the entangled network is estimated as the elastic energy of all its stretched entangled segments, A more rigorous estimate for the free energy of the non-affine tube model with the number of monomers of the chain N N e gives [29] F e V = ρkT This expression includes an additional to Equation (48) contribution ∼1/λ α describing the compression of the chain in the tube.

Slip-Tube Model
This model generalizes the nonaffine tube model by taking into account the slippage of chains along the tube contour. In the case of anisotropic deformation of the network, its chains are more stretched in the direction of maximum extension relative to other directions. The increased chain tension in this direction draws the chain from the tube segments along these directions. The redistribution of the chain length between different sections of the tube is taken into account by renormalization of parameters of the non-affine tube model for each direction of deformation The parameters g α are normalized by the condition of preserving the full length of the chain, and in the case of uniaxial deformation, Equation (3), there is only one independent redistribution parameter g z : In addition to the above renormalization of Equation (49), the free energy of this model has an entropic contribution taking into account the chain slippage along the tube, The parameters g α are found by minimizing the free energy (62), and the minimization equation reads as The elastic stress of the slip-tube model is The effect of entanglements substantially depends on the degree of network stretching, and the chains of strongly stretched networks become effectively phantom. Such a crossover from entangled to affine behavior is usually described using the assumption of additivity of the phantom and affine contributions to the free energy: This assumption also leads to the additivity of the corresponding contributions to the elastic tensor, In the case of uniaxial network deformation, Equation (3), a convenient representation of stress is the Mooney ratio of the total stress σ zz − σ xx to the functional dependence λ 2 − λ −1 predicted by phantom network models: As was shown in Ref. [3], the Mooney plot of the experimental data is well described by the slip tube model.

Crossover from Unentangled to Entangled Regimes
The Mooney-Rivlin dependence (58) is usually fitted by the phenomenological equation It is usually assumed that the parameter 2C 1 is only due to chemical cross-links, while the parameter 2C 2 includes all of the entanglement contributions. According to the non-affine tube model, part of the entanglement contribution is also included in the parameter 2C 1 [6].
In this section, we study the concentration dependence of the Mooney-Rivlin coefficients. Although the free energy additivity approximation (56) can be used to describe the polymer network in the melt, a refined theory should be developed to find the Mooney dependence in the entire concentration range. The importance of such a study is emphasized by the fact that the networks prepared near the crossover from unentangled to entangled regimes are a fairly typical case.
Consider a network, consisting of ν Gaussian chains with N monomers whose endpoints are displaced affinely with the global deformation of the network. We assume that each chain of this network is confined in an effective "non-affinely deformed" tube, described by the virtual chain's potential (43). The free energy of this model is calculated as the sum of contributions of all its strands and is found in Appendix C: where k α = N/ (N e λ α ). In the "phantom" limit N N a f f α = λ α N e Equation (60) turns into the free energy of affine network model, Equation (4). In the limit of strongly entangled networks, N N a f f α , this expression reproduces the free energy of the non-affine tube model, Equation (49). The first term in Equation (60) can be rewritten as where R 0 = bN 1/2 is the chain size under network preparation conditions and R f l,α is the fluctuation size in the α-direction, equal to R 0 in the phantom regime and the tube diameter a α in the entangled regime. The logarithmic term in this expression describes the entropy of the chain compression in the tube.
By renormalizing the parameters of this model (50) to take into account the effect of chain redistribution in the entangled tube, and adding the contribution of entropy of slippage along the tube (53), we arrive at the expression for free energy Solving equations (54) for this model, we find the stress (55) in the polymer network, which determines the Mooney-Rivlin dependence (58). The parameter C 2 (φ) of this dependence is determined by the expression We plot the dependences C 1 (φ) and C 2 (φ) (63) for λ −1 = 0.7 and N/N e = 2.4 in Figure 4. The parameter C 1 of the Mooney-Rivlin dependence (59) weakly depends on the polymer volume fraction φ. In accordance with the experimental data [30], the dependence of the parameter C 2 on φ is almost linear, and C 2 vanishes at φ 0.2. The vanishing of C 2 at finite φ is related to the transition between phantom and entangled regimes. When the polymer swells, its chains move away from each other, which leads to a weakening of the effective topological potential that holds the chains in the tube. Chain fluctuations in this potential increase with the swelling and can reach the fluctuations of the phantom chain. This concentration corresponds to the disappearance of the contribution of entanglements to the network elasticity and the vanishing of the Mooney-Rivlin coefficient, C 2 = 0. benzene carbon tetrachloride petroleum ether (60/80) nitrobenzene decane Figure 4. Dependence of the Mooney-Rivlin coefficients C 1 (φ) (upper curve) and C 2 (φ) (lower curve) on the volume fraction of polymer φ: theory (solid curves) and experimental data (points, Kg/cm 2 ) for swollen rubber [30].

Heterogeneities in Polymer Networks
Network inhomogeneity is a common feature of polymer networks and gels due to the randomness of the crosslinking process and the presence of additional topological defects such as dangling chain ends, cross-linker shortcuts, and chains forming loops. The origin of nanostructural inhomogeneities and their characterization by light, neutron, and X-ray scattering as well as by NMR spectroscopy and optical, electron, and X-ray microscopies is reviewed in Ref. [31], and the main methods of their study are outlined in the review [32]. The first attempts to take such heterogeneities into account were made by using the model of randomly cross-linked networks containing fractal regions, such as percolation clusters [33]. Networks with fractal heterogeneities can be swollen and deformed by unfolding the fractal regions without significant elastic entropy penalty [34]. It is shown that strong heterogeneities lower the shear modulus of the network if a part of strands is so short to be considered rigid and not able to deform. Networks with a large number of such very short strands have higher breaking energies.
The effect of entanglements on heterogeneities in polymer networks is studied in Ref. [35] using molecular dynamics simulations of polymer networks made by either end-linking or randomly crosslinking a melt of linear precursor chains. The end-linking leads to nearly ideal monodisperse networks, while random cross-linking produces strongly polydisperse networks. It is shown that the microscopic strain response, the diameter of the entanglement tube, and stress-strain relation weakly depend on the linking process by which the networks were made.
The replica method was used in Ref. [28] to calculate the characteristic size and amplitude of the spatial nonuniformities of the network due to defects of its structure and topological restrictions. Using this method it is shown that inhomogeneities can arise as consequences of a stretching of polymer networks [36]. Although the replica theory [8] provides a complete solution of the statistical mechanics of polymer gels, it uses replica trick which is unfamiliar to the majority of people in the polymer community. A more intuitive phenomenological approach capturing all the main physical ingredients of the complete theory is developed in Ref. [37].

Theory of Heterogeneities in Polymer Networks
In this section, we review the main results of the theory of random spatial heterogeneities developed in swollen and deformed polymer networks [38]. The non-triviality of this problem stems from the fact that information about network structure is "encrypted" in the pattern of cross-links joining polymer chains, which represent a very small fraction of the network volume. The initial crosslinking pattern is reproduced only partially due to thermal fluctuations arising in the new equilibrium state after crosslinking. Conformations of polymer strands in such networks with fixed topological structure can be varied in a wide range depending on experimantal conditions. The density profile of monomers in the polymer network can be recovered from the Fourier component of the deviation of the density from its average value ρ: Here, δρ(x) are random deviations (due to thermal fluctuations) of the density from the equilibrium density profile ρ eq (x) describing spatial inhomogeneities in polymer networks. The fluctuation contribution to the free energy of the network with a given distribution of monomer density can be represented as the sum of the entropy contribution and the (osmotic) contribution of volume interactions ∼ v: Here n (q) is the density profile in the "elastic reference state" [37], maximizing the entropy of the polymer network. The density n (q) is determined by the molecular structure of the network and it vanishes in the short-wavelength limit q R −1 , since solid-state degrees of freedom are determined only at length scales exceeding the monomer fluctuation radius R. On a smaller scale, the gel has liquid degrees of freedom, contributing to the temperature structural factor The first term in square brackets determines the contribution of the liquid degrees of freedom of the polymer network. The dimensionless wavevector Q = Rq is normalized by the monomer fluctuating radius, and the vector λ · Q has components λ x Q x , λ y Q y , λ z Q z . At large Q 1 the term Q 2 /2 in the denominator of the first term of the right hand side of Equation (66) gives the usual Lifshitz entropy of polymer solutions,g (q) = 2ρN/Q 2 [39]. The term (4Q 2 ) −1 , first introduced by de Gennes for heteropolymer networks [40], describes the suppression of density fluctuations on length scales larger than the monomer fluctuation radius R. The second term in square brackets in Equation (66) determines the contribution of solid-state degrees of freedom of the polymer network. It remains finite in the long-wavelength limit and retains its angular dependence on the anisotropic deformation even in the limit q → 0.
Equilibrium monomer density profile is found by minimizing the free energy (65) As can be seen from Equation (67), in good solvent with positive excluded volume v > 0, the equilibrium density profile is more homogeneous than that of the corresponding elastic reference state. The excluded volume interaction suppresses not only thermal fluctuations in the polymer network, but also frozen-in spatial inhomogeneities. Although the static spatial fluctuations are permanent inhomogeneities, they can reversibly increase and diverge at the spinodal line, at which vg (0) = −1. This observation is experimentally confirmed in Ref. [41].
In (ergodic) heterogeneous systems, there are two types of averages, i.e., the thermal or time averages and ensemble or space averages, denoted by X and X, respectively. The Fourier component of the thermal correlator of density fluctuations δρ (x) = ρ (x) − ρ eq (x) is found by averaging with the Gibbs weight e −F [ρ]/kT , Equation (65): Therefore, in addition to permanent spatial heterogeneities, a polymer network undergoes thermal dynamical density fluctuations which diverge at the spinodal line.
All information about the heterogeneity of the molecular structure of the network is contained in the monomer density n (q) in the "elastic reference state", which can be considered as a random variable characterized by a correlator Therefore, it is important to understand the physical meaning of the different terms of this expression: Each of two factors (1 + Q 2 ) −1 describes the thermal "smearing" of the density response to random stresses in the polymer network. The first term in square brackets comes from the correlator of frozen random stresses, which exist even in a homogeneous polymer network. Such stresses appear due to heterogeneities in the distribution of cross-links. The second term in brackets in Equation (69), which is proportional to the correlator of the affine deformed density pattern in the reference state S (0) (λ · q), strongly depends on network preparation conditions.
(1) In the case of cross-linking in a melt, density fluctuations are suppressed andS (0) (q) = 0. The fact that ν (q) does not vanish upon swelling of such a network, despite the absence of density fluctuations at the preparation conditions, means that there are still inhomogeneities in the crosslink density, which manifest themselves upon swelling.
(2) In the case of instant cross-linking of a semi-dilute polymer solution,S (0) (q) is given by the Ornstein-Zernicke expressioñ which is finite for any second virial coefficient v 0 at network preparation conditions.
(3) In case of equilibtium chemical cross-linking, the total structure factor of the gel in the state of preparation is determined by the expression for the polymer liquid, Equation (70), with renormalized The decrease in v (0) is due to the fact that the monomers forming cross-links give an additional attractive contribution to the effective ren . The amplitude of heterogeneities increases when approaching the cross-link saturation threshold at v (0) ρ (0) N = 1 at which both the structure factor S (0) (q → 0) and the characteristic size of heterogeneities at preparation conditions The theory also predicts the appearance of a maximum degree of spatial gel inhomogeneity at a critical polymer network concentration, as confirmed by experiments on PAAm gels [42,43].

Scattering Intensity
The scattering intensity on wavevector q is proportional to the structure factor, which is given by the sumS of contributions of the thermal fluctuations, Equation (68), and the inhomogeneous equilibrium density variations due to defects of the topological structure of the network, where bar means ensemble average. The experimental data can be visualized using contour plots of neutron scattering from random inhomogeneities of network structure in anisotropic-deformed swollen gels [44]. Under uniaxial gel stretching, an increase in the scattered intensity on small wave vectors q in the stretching direction is observed, which is enveloped by elliptical patterns at larger q values with a maximum oriented normal to this axis. This behavior is opposite of what is expected in theories assuming only thermal fluctuations and is called "abnormal butterfly patterns". The elliptical patterns at large wavevector q originate from the correlator of static inhomogeneities, Equations (72) and (69), which contains the term S 0 (λ · q) that "remembers" the affinely deformed inhomogeneous structure of the network. The butterfly patterns along the stretching axis in the small q range is due to strong angular dependence of the thermal structure factor, (68). This function comes into numerator of the thermal correlator in Equation (68), describing "normal butterfly patterns". At high cross-link concentrations the correlator of static inhomogeneities gives the main contribution to the scattering intensity. Since the functiong (q) is included in the denominator of Equation (72), this expression describes the so named "abnormal butterfly patterns". As a result, under uniaxial extension, a crossover from normal to abnormal butterfly scattering patterns occurs with an increase of the strength of inhomogeneity or the swelling ratio [45]. In addition to such butterfly scattering patterns, the theory also predicted "Lozenge" patterns if only part of all network chains, i.e., their deuterated fragments, could scatter neutrons, in accordance with scattering experiments [46].
The static inhomogeneity in poly (N-isopropyl acrylamide) gel (PNIPA) has been investigated in Ref. [47] by the methods of small-angle neutron scattering (SANS) and neutron spin echo. The obtained SANS scattering amplitudeS (q) was successfully decomposed into the thermal and static components, respectively,D (q) andC (q) in Equation (71). It was revealed thatC (q) becomes dominant in the q-region, where the abnormal butterfly scattering under stretching is observed. As the temperature increases toward the temperature for volume phase transition,C (q) approximated by the square of the Lorentzian shape increases more drastically thanD (q) of the Lorentzian shape. These experimental findings are also well described in the theoretical framework of this section.
In Ref. [48] gels prepared by two cross-linking methods were studied using SANS technique. One is chemical cross-linking with BIS and the other is gamma-ray cross-linking of a PNIPA solution. It is shown that the degree of the inhomogeneity is much larger in chemically cross-linked gels than in the gamma-ray cross-linked gels. The experimental data are in quantitative agreement with predictions of the theory of scattering intensityS (q) on chemically and instantaneously cross-linked gels, respectively.

Heterogeneities in Charged Gels
A study of the structure factor of weakly charged polyelectrolyte gels under uniaxial stretching was performed by Mendes et al. [49], who observed after introducing ions to the gel, the disappearance of the "butterfly pattern" in the small angle scattering intensity, as well as an increase in the scattering intensity in the direction perpendicular to the gel stretching. The origin of this maximum has been elucidated in SANS experiment by Shibayama et al. who studied the deformed state of weakly charged polymer gels PNIPA/AAc [50]. In this experiment, an anisotropic scattering maximum is observed, which indicates that the spatial distribution of the charged groups changes as a result of gel deformation and therefore is strongly coupled with the static inhomogeneities.
All the observed patterns of SANS intensity were well reproduced using a generalization of the above theory to the case of charged polymer networks [51]. The only effect of electrostatic interactions is to replace the excluded volume parameters v (0) and v with effective virial coefficients for both the final state and for the state of preparation: is the inverse of the Debye-Hückel structure factor, κ −1 is the Debye screening length, α is the degree of ionization (fraction of charged monomers) and l B is the Bjerrum length. s (0) DH (q) is obtained by substituting κ = κ (0) and α = α (0) into the above equation. As shown in Ref. [52], the above theoretical prediction for S (q) well reproduces the observed scattering intensity functions of wekly charged PNIPA/AAc gels.
In general, an introduction of cross-links to a polymer solution leads to an increase in the scattering intensity due to static inhomogeneities. However, a reverse phenomenon, called the "inflection" in scattering intensity, was predicted by the above theory [53] and observed in weakly charged gels and polymer solutions [54]. While the gel becomes more inhomogeneous with increasing the degree of cross-linking in a good solvent, the inhomogeneities can be suppressed in a poor solvent, although in a relatively small regions of cross-link concentration. However, this phenomenon is interesting due to its physical significance.

Good Solvent: Scaling Approach
Earlier in this paper, we considered the networks with Gaussian chains obtained in the melt or the θ-solution of linear chains. In general, polymer networks can also be obtained by crosslinking semi-diluted polymer solutions. The networks can be placed in a good solvent, in which the polymer chains swell compared to the case of theta solvent. This section explains the basic ideas of the scaling approach to describing networks prepared/swollen in a good solvent.
The mean-field approach can be adapted to descripbe the elasticity of such gels [55] using the well-known de Gennes blob picture of semi-dilute solutions [56]. The key idea is the spatial scale separation: while static density inhomogeneities exist only on scales comparable to or larger than the monomer fluctuation radius R, thermal density fluctuations are dominated by smaller scales and are quite similar to those in semi-dilute polymer solutions.
Consider a gel prepared by random cross-linking of chains in a semi-dilute polymer solution in a good solvent at the monomer density ρ (0) that is swollen to density ρ < ρ (0) . The monomer density fluctuations should be taken into account both in the initial (where the gel was prepared) and in the final (where it is being studied) states of the gel. On length scales shorter than the correlation lengths in these states, density fluctuations are large, and the gel behaves like a polymer solution ("liquid-like" regime). On scales exceeding the blob sizes (the corresponding wave vectors q (0) = 1/ξ (0) and q = 1/ξ), density fluctuations are small, and the mean-field description with appropriately renormalized parameters can be used [8] b ren /b ren where the subscript ren stands for renormalized values that differ from the "bare" ones.
The renormalized second virial coefficients v (0) ren and v ren are: Equations (73) and (74) complete the renormalization of the mean-field theory: in order to describe a gel in a good solvent on length scales larger than the thermal correlation length, one only have to replace the bare parameters in the previously derived expressions for the free energy, correlation functions, etc., by their renormalized values.

Amplification of Cross-Linking Density Pattern
Is it possible to "write" some "useful" information in a polymer network cross-link pattern, and under what conditions it can be "read" back? The answer to this question is given in Ref. [57]. After the formation of a homogeneous (on length scales large compared to its "mesh" size R) network, large-scale patterns can be generated in it by further cross-linking followed by swelling (and possibly stretching) of the network, which leads to a nonuniformly swollen gel. This additional cross-linking can be done by adding light-sensitive cross-links to the transparent network. By focusing the laser beam in the regions inside the gel, it is possible to "write" information into the gel structure in the form of 2D or 3D patterns of cross-link density.
It was shown that although such information is hidden at the preparation conditions, it can be reversibly "developed" by gel swelling, since unobservable variations of the cross-link density in the melt are transformed into observable variations of monomer density in the swollen gel. The gel regions with increased cross-link concentration can be considered as inclusions with enhanced elastic modulus. It has been shown that in swollen gels that stretch isotropically upon absorption of the solvent, the observed monomer density pattern is an affinity stretched initial cross-linking pattern. Such gels can serve as a magnifying glass that enlarges the initially written pattern without distorting its shape. The corresponding magnification factor can be very large in the case of super-elastic networks.
A strong enhancement of the contrast between the high and the low monomer density regions can be obtained by placing the gel in a poor solvent with a negative second virial coefficient v < 0. The observable image becomes distorted, especially near the edges and corners of the pattern. Gel boundaries should be fixed in order to avoid its contraction. Gel contraction can also be prevented by focusing a laser beam only on a part of the localized pattern and heating it, resulting in a local change of the quality of solvent.

Discussion
Despite its more than half-century history, even now the theory of polymer network elasticity asks more questions than it gives answers. In this work, we tried to bridge the gap between the main developed approaches, which allowed us to give answers to some of these questions.
In real networks, there are both topological defects and typical loops, which are not sparsely distributed. Using the replica method, the cumulative effect of a large number of strongly overlapping typical loops of finite size can be described in the effective mean-field approximation. This approximation works for polymer networks due to large overlap parameter of network strands, P (0) 1, see Equation (12). In the replica approach, the properties of polymer networks are described by the liquid-solid (sol-gel) order parameter ϕ 1 (ς), which is actually not a parameter, but a function of the variable ς, and is determined by a complex integro-differential equation. The only, but very important exception is the elastic modulus, which is expressed through the value of the order parameter at ς = 0, and the corresponding equations for ϕ 1 (0) become algebraic. This greatly simplifies the calculations of the elastic modulus of the network with an arbitrary number of cyclic fragments of finite size. In contrast to the classical theory of phantom networks, the mean field of loops, ϕ 1 (0) explicitly depends on the excluded volume parameter v (0) and the density of monomers ρ (0) at network preparation conditions. This dependence takes into account the limitations imposed by the packing of highly overlapping typical loops in real 3D space on the molecular structure of the network being formed.
The resulting elastic modulus of real networks has two main contributions, see Equation (37). The contribution of elastically effective network strands is taken into account in the classical model of phantom networks. The classical expression for this contribution is renormalized due to the presence of topological defects in the network-primary loops and cyclic fragments of finite size. The second contribution to the elastic modulus is always negative and gives an impact of a large number of typical loops of the network. It describes the effect of shunting of the elastically effective chains by finite size polymer loops and depends on the interaction of monomers at the preparation conditions of the network.
Both contributions to the network elastic modulus can be represented by the generalized combined chain model with an additional effective chain (see Figure 2b) that accounts the cumulative effect of finite size loops in real networks. The virtual chains transmit local stresses from the network fragments to the solid-state degrees of freedom of such soft solids. We established the connection of the distribution function of the lengths of virtual chains with the order parameter of the replica network model and calculated this distribution function for the exponential distribution of the lengths of the network strands.
We also proposed a generalization of the slip tube model of entangled polymer networks, which allows us to describe the crossover between the entangled and phantom regimes of network swelling. The dependence of the Mooney-Rivlin parameters C 1 and C 2 on the polymer concentration calculated for this model is in agreement with the experiments.
Although from the very beginning the heterogeneities were recognized as one of the most essential features of gels [58], it took long time to formulate their theoretical description due to complexity of this problem. Up to today, the understanding of the gel structure has greatly improved owing to both theoretical developments and a large number of experimental studies. Effect of cross-links on heterogeneous structure of polymer gels, abnormal butterfly patterns, microphase separation, and so on, are well understood with the aid of the theretical approaches discussed in Section 4.
The list of remaining questions is much longer: Polymer networks can be obtained by crosslinking semidiluted polymer solutions. How does a strong nonlinear deformation of such networks occur in the phantom regime? How do such networks fracture when they are strongly stretched? The molecular theory (Lake-Thomas model) of the fracture of dry networks is constructed in Ref. [59]. What happens when compressing such networks? What is deformation mechanism of heterogeneous super-tough networks, kinetics of crack growth and phase transitions in such networks? Stress-induced microphase separation in multicomponent networks was studied in [60]. What is the effect of microphase separation on the elasticity of such networks?
These and many other questions are still awaiting answers.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Replica Model of the Network
Consider a network prepared by random end-linking a solution of polydisperse linear chains by f -functional cross-links. Following the idea of ϕ 4 formulation of the excluded volume problem [56], we introduce n → 0 component field, ϕ (x) in the replica space with components ϕ i (i = 1, . . . , n), which is related to the monomer density in this space as Equations (25) and (A6)-(A8) completely determine the distribution function p (n) of the number of monomers n of virtual chains for a given distribution function P (N) of the number of monomers N of real network strands. The only problem is solving this system of nonlinear integral equations. Here the replica method comes to the rescue, see Appendix A.

Appendix C. Free Energy of Entangled Networks
The network free energy for given positions of attachment points X(s) of virtual chains to the affinely deformed background is The simplest way to calculate the above integral is to use the normal mode expansion of the chain trajectories, which satisfy the boundary conditions x α (N) = λ α X α (N) and x α (0) = λ α X α (0) (X (s) are the positions of attachment points of virtual chains in the undeformed network): The mode with ω = 0 is excluded from the summation because it corresponds to the translational displacement of the chain, which is prohibited in the network. Substituting the above expressions into Equation (A10) we get Changing the integration variables,x i (ω) =ũ (ω) + iṽ (ω), we find Here sums are taken only over positive ω since the only positive frequency componentsũ (ω) andṽ (ω) can be taken as independent variables: the condition for the reality of the function x(s), Equation (A12), imposes restrictionsũ (−ω) =ũ (ω) andṽ (−ω) = −ṽ (ω) on the real and imaginary parts ofx (ω). This free energy should be averaged over positions of attachment points, X(s), with correlation functionsX α (ω)X β (−ω) = Averaging the free energy (A14) with the correlation function, (A15) and changing the integration variablesX(ω) → (R,p (ω) ,q (ω)) when integrating over X (ω) we find Here all the terms which do not depend on the extension coefficients λ α are absorbed in Const. The sum and the product can be expressed in terms of hyperbolic functions, see Equation (60).