Coarse-Grained Lattice Modeling and Monte Carlo Simulations of Stress Relaxation in Strain-Induced Crystallization of Rubbers

Two-dimensional triangulated surface models for membranes and their three-dimensional (3D) extensions are proposed and studied to understand the strain-induced crystallization (SIC) of rubbers. It is well known that SIC is an origin of stress relaxation, which appears as a plateau in the intermediate strain region of stress–strain curves. However, this SIC is very hard to implement in models because SIC is directly connected to a solid state, which is mechanically very different from the amorphous state. In this paper, we show that the crystalline state can be quite simply implemented in the Gaussian elastic bond model, which is a straightforward extension of the Gaussian chain model for polymers, by replacing bonds with rigid bodies or eliminating bonds. We find that the results of Monte Carlo simulations for stress–strain curves are in good agreement with the reported experimental data of large strains of up to 1200%. This approach allows us to intuitively understand the stress relaxation caused by SIC.


Introduction
Natural rubbers undergo strain-induced crystallization (SIC), which attracts a lot of attention and has been studied extensively [1][2][3][4]. When the rubbers are stretched by an external tensile force, the entangled polymers partly change into an aligned state, and, then, this aligned state gradually converts into a crystalline state (Figure 1a,b). Moreover, if the stress is released after stretching or loading, one can also observe that the crystallized part immediately starts to melt. In this recovery or unloading process, hysteresis can be observed in the stress-strain curve and crystallization ratio (Figure 1c). This hysteresis indicates that the stress decreases due to the so-called crystallization-induced strain relaxation [3], which can be called stress relaxation. Such a stress relaxation is considered to be an equilibrium property because the stress-strain curve during the unloading process is close to the equilibrium curve. However, the stress relaxation is contradictory to our intuitive understanding that the mechanical strength is expected to be increased by the crystalline part of polymers, of which the mechanical strength is very large compared to the amorphous part, at least for intermediate strain region. Therefore, the stress relaxation by SIC should be studied in more detail to clarify how it can be understood by a simple modeling technique on the basis of statistical mechanics. Ref. [5] τ χ χ loading unloading A B Figure 1. (a) The h state of rubber, where the polymer directions are isotropic, (b) extended or elongated state of rubber, where the polymer directions align along the tensile force direction and strain-induced crystallization (SIC) starts to appear, and (c) an illustration of stress-strain curves and crystallization ratio χ for loading and unloading processes influenced by SIC. Small circles in (a,b) represent cross-linkers. In (c), the letter A (B) denotes the point where crystallization starts (terminates) in the loading (unloading) process. The stress τ(MPa) and the crystallization ratio χ(%) in (c) are taken from experimental data in [5].
It is also interesting to see that the changes in state observed in the SIC process share the universal property with the first-order transition between gas and liquid; for example, [3]. From this first-order nature of SIC, we understand that the plateau behavior in the stress-strain curve for the unloading process has the same origin as that of the plateau in the pressure-density curve of the gas-liquid transition. On the other hand, the difference between SIC and the gas-liquid transition is that the crystallization is seen only partly in the case of rubbers, and this partial crystallization seems to occur due to the fact that the entangled polymers cannot be disentangled by the simple loading process from the topological constraint.
Therefore, in general, the effect of SIC on the stress-strain curves is very difficult to implement in the model [6]. In fact, it is easy to understand that the amorphous state is mechanically very different from the crystalline state [7][8][9][10][11][12]. Since the stress and strain of the crystalline state are expected to be completely different from those of the amorphous state, the mechanical response to external force is also completely different between the two states. For this reason, the tensile energy or Hamiltonian of the amorphous state cannot be shared by the crystalline state. This is one of the reasons why SIC is hard to study from the viewpoint of statistical mechanical modeling, where the Hamiltonian for each microscopic state should be defined explicitly.
In this paper, we use the tensile energy of the amorphous state of polymers simply by discarding the implementation of energy for the crystalline state. The tensile energy of the amorphous state is the so-called Gaussian bond potential, of which the one-dimensional version was originally assumed in the Gaussian chain model for polymers, where "bond" corresponds to the segment of the linear chain [13]. The two-dimensional (2D) version of the Gaussian bond potential is used for the triangulated surface models of membranes [14][15][16][17], and the 3D extension is straightforward with a tetrahedral lattice. In this paper, 2D and 3D versions of the Gaussian bond potential are assumed for our new SIC modeling. For these potential energies, the crystalline state can also be implemented simply by replacing elastic bonds with rigid bonds or empty bonds, where the rigid bond is a three-dimensional rigid body of fixed length and has no tensile energy, and the empty bond also has no tensile energy; these models are called the "rigid bond model" and "empty bond model", respectively, in this paper. These two models are close to each other in the sense that no tensile energy is defined on the bond corresponding to the crystalline state.
It is well known from the thermodynamics of rubbers that rubber elasticity comes mainly from entropy [18,19]. Since this entropy contribution is a major part of the elasticity, the rubber elasticity is called entropy elasticity. However, we study rubber elasticity by statistical mechanical models in which two different microscopic variables are assumed, and we will not go into the details of thermodynamic theory.

2D Lattice
A 2D triangulated lattice was originally constructed for the surface model of membranes [14][15][16][17]. Here, we use cylindrical lattices for the calculation of surface tension by fixing the boundary of the cylinder, which will be shown in the following subsection. The lattice size N, which is the total number of vertices including the boundary vertices, is given by N = L 1 (L 2 −1), where L 1 and L 2 are the total number of vertices on the edges of the rectangular plate ( Figure 2a). The height H 0 and the diameter D of the cylinder are assumed to be the same such that where a is the lattice spacing or can be called the edge (or bond) length of the initial lattice. The lattice spacing is fixed to a = 1 and eliminated from the expressions of bond length henceforth; however, this spacing is restored if the units of surface tension are changed from the simulation units to the physical units in a later section. The diameter D is given by πD = (L 2 −1)a. Therefore, from D = H 0 , the relation L 2 = ( √ 3π/2)(L 1 −1)+1 is the condition for the triangles to be the regular triangle. A snapshot of cylinders of size L 2 = 14, L 1 = 6, and N = 78 is shown in Figure 2b, where the condition L 2 = ( √ 3π/2)(L 1 −1)+1 is not always exactly satisfied because the right-hand side is not an integer.

;ĂͿ ;ďͿ
(a) A cylindrical surface is made of a rectangular surface of size L 1 ×L 2 , where L 1 and L 2 are integers, and (b) a triangulated cylinder of size N = 78, where the total number of boundary vertices is N bnd = 2(L 2 −1) = 26. Small numbers L 1 (= 6) and L 2 (= 14) are assumed here to visualize triangles clearly. The height H 0 and the diameter D are given by Equation (1) with the edge length a.

Rigid Bond Model and Empty Bond Model
Let r i (∈ R 3 ), (i = 1, · · · , N) be the vertex position of the lattice (Figure 3a), ij = | r i − r j | be the bond length (Figure 3b), and n i be the unit normal vector of the triangle (Figure 3c). These variables r i are integrated into the partition function Z for 2D models (here, we discuss 2D surface models), and Z is given by where the symbol 0 in Z( 0 ) indicates that the model implicitly depends on a constant 0 , which will be used to define the crystalline bond below. The symbol d r i denotes the 3D integration for the internal vertices. The other N bnd i=1 d r i denotes 2D integration for the boundary vertices on the upper and lower boundaries with a periodic boundary condition such that where Z U i and Z L i denote the Z component of the vertex position r i = (X i , Y i , Z i ) on the upper and lower boundaries (no confusion is expected for the symbol of the partition function and that of the coordinate symbol). The second equation in Equation (3) denotes a constraint that the diameter of the boundary circle is given by D (Figure 2a). The reason why N bnd i=1 d r i becomes a 2D integration is that the boundary vertices are allowed to move not only along the boundary circle of fixed diameter D, but also along the Z direction due to the periodic boundary condition. For the 3D thick cylinder model, N bnd i=1 d r i becomes a one-dimensional integration because no periodic boundary condition is imposed. The symbol σ in Equation (2), which has values in {0, 1}, is another dynamical variable defined on the bond i. The value of σ represents whether the bond i is rigid or amorphous: The rigid bond is a rigid body, which is allowed to move by 3D rotation and translation. In other words, two vertices i and j, which are connected by a rigid bond, move with a constraint ij (= | r i − r j |) = c ij , where c ij is described in the following paragraph. In contrast, this constraint on ij is not imposed on any two vertices connected by an amorphous bond. Another model is the "empty bond" model, which is defined by replacing "rigid" with "empty" in Equation (4) such that , (empty bond model).
The variable σ i in Equations (4) and (5) can be changed from σ i = 0 to σ i = 1 when the bond length i is i > 0 using the constant 0 in Z( 0 ) according to the Monte Carlo update procedure described in a later section. If σ i is changed from σ i = 0 to σ i = 1, then the bond length i is fixed to the length c i , which is the length of bond i just before the change in σ i in the rigid model. In the case of the empty bond model, the length i of the crystalline bond is changeable to a new length i if i < c i , and this c i is also changeable depending on the update of σ i , which is the same as in the case of the rigid bond model. The constraint c i > 0 is always satisfied in both models. The empty bond is almost close to the rigid bond except that the length i is changeable to i < c i . This fact implies that the empty bond length i is also prohibited from extending to i > c i . The bond configurations for amorphous, rigid, and empty bonds are illustrated in Figure 4a-c. We should comment that the total number of degrees of freedom of variable r is slightly reduced by crystallization. Indeed, in the rigid bond model, a crystalline bond is allowed to move by three-dimensional translation and two-dimensional rotation, where the rotational degrees of freedom are given as 2 because the bond is a one-dimensional object. Therefore, the 6 degrees of freedom, which is the sum of the degrees of freedom for translations of the terminal vertices of the bond, is reduced to 5 by crystallization. Thus, the total number of degrees of freedom is reduced by N cr , which is the total number of crystalline bonds, in the rigid bond model. In the case of the empty bond model, the total number of reductions in the degrees of freedom of variable r is effectively considered to be N cr /2 because the crystalline bond length is allowed to shrink even though the extension is not allowed, and, for this reason, the reduction per crystalline bond is roughly estimated to be 1/2.
The Hamiltonian for the rigid and empty bond models is given by and where the first two terms S 1 +bS 2 and the third term κS 3 correspond to the tensile energy and the bending energy, respectively, and the final two terms U 1 +U 2 are constraint potentials. The symbol σ in the definition of ij is the variable σ on the bond connecting the vertices i and j. No difference between the rigid and empty bond models is in the definition of ij except that the length ij of crystalline bond i j is allowed to be ij < c ij in the case of the empty bond model. In contrast, the length of bond i j is prohibited from being ij > c ij for σ ij = 1 in both models. In this sense, this definition of ij is specific to the models for SIC in this paper. In the bending energy S 3∆,i (i = 1, 2, 3), σ is defined on the bond shared by triangles 0 and i (see Figure 3c). This definition of the bending energy is also specific to the models for SIC in this paper. These specific definitions of ij and 1− n i · n j come from the fact that the crystalline state is implemented in the context of the Gaussian bond model for polymers [13].
Due to the factor 1/2, the Gaussian bond potential S 1 is identical to S 1 = bond i 2 i , which is given by the sum of bond length squares 2 i . The reason why the sum ∆ is used for S 1 instead of bond i is because the discrete form of the quadratic term S 2 is naturally given by ∆ , and the same summation convention is used for S 1 . The numerical factor 1/2 in S 2 is determined such that it is identical to the factor in S 1 .
Note that U 1 is the constraint potential for i in both models. From this constraint U 1 , the length i of the empty bond is only allowed to be i ≤ c i , as mentioned above; however, this is strictly prohibited for the rigid bond because i should always satisfy i = c i from the definition of i . The second constraint U 2 prohibits the crystalline bond i from being consecutively connected (Figure 5a), or, in other words, every crystalline bond i must be isolated or linked to the amorphous bond j (Figure 5b). This constraint U 2 also plays a role in the constraint on the total number of crystalline bonds because the maximum number of crystalline bonds is limited by a lattice geometric (or topological) reason.
As a consequence of the constraint U 2 , no constraint is imposed on the upper limit for the crystallization ratio χ such that χ ≤ χ max , where χ is defined by where i σ i is the total number of crystalline bonds and N B (= i 1) is the total number of bonds. The value of χ max depends on the dimension of the lattice, that is, whether it is a 2D or 3D lattice; however, in both cases, χ max is reasonable in the sense that it is close to the experimental value, such as χ max (Exp) 0.1 [1][2][3][4], which will be shown in the presentation section. The mean value Q of physical quantity Q( r, σ) (where the symbol Q denotes an arbitrary quantity) is given by We simply write the mean value Q by removing the symbol · henceforth. To calculate these multiple integrations and the sum σ over all possible configurations of σ, we use Metropolis Monte Carlo (MC) simulations, which are described in the following section. Note that all physical quantities are calculated with the dynamical variables r and σ, which are the vertex position r = { r 1 , · · · , r N } and the spin variable at the vertices σ = {σ 1 , · · · , σ N }, respectively.
;ĂͿ ;ďͿ ሺൈሻ ሺ ሻ Figure 5. (a) A configuration that is not allowed because the two bonds of σ = 1 are connected, and (b) a configuration that is allowed because the two bonds of σ = 1 are separated. The maximal value χ max of the crystallization ratio of Equation (8) is not given as an input parameter in the simulations, but it is automatically determined by the constraint illustrated in the figure.
The 2D continuous forms of the Hamiltonian corresponding to the discrete S 1 and S 2 in Equation (6) are given by where g ab is the inverse of a metric tensor g ab given by the 2 × 2 matrix, and g is the determinant.
To obtain the discrete Hamiltonian, we simply assume that g ab = δ ab , which is called the Euclidean metric. By the replacements ∂ 1 r(= ∂ r/∂x 1 ) → r 2 − r 1 and ∂ 2 r → r 3 − r 1 (see Figure 3b), and by the symmetrization of the obtained expression of the Hamiltonian due to the three possible local coordinate origins in a triangle, we have the discrete S 1 and S 2 in Equation (6). More detailed information on the discretization is given in Ref. [20].

3D Model
The three-dimensional (3D) continuous Hamiltonian is obtained from the 2D version in Equation (10) simply by replacing the 2D integration √ gd 2 x with the 3D integration √ gd 3 x.
In the 3D case, the metric tensor g ab is also replaced by a 3 × 3 matrix. Thus, these 3D continuous S 1 and S 2 can also be converted to discrete forms on the 3D lattice, which is a thick cylinder (Figure 6a) discretized by tetrahedrons ( Figure 6b). On this thick cylinder, the vertices of tetrahedrons are distributed only on the inner and outer surfaces, and, therefore, the thickness is negligible compared with the diameter. Note also that in the 3D case, only the empty bond model is studied. The reason is that, as we will see in the presentation section, there is almost no difference between the results of the 2D empty and rigid models, or the 2D empty model is slightly better for the stress relaxation behavior, and, moreover, the empty model is more simple than the rigid model in its definition. The 3D version of the constraint potentials U 1 and U 2 has the same expressions as those in Equation (7), and, hence, their expressions are not written below.
where tet in S 1 and S 2 denotes the sum over the tetrahedrons (Figure 6a). The symbolN is the mean value of the total number of tetrahedrons around a bond and is given bȳ where tet(i) is the total number of tetrahedrons around bond i, and N B is the total number of bonds. Due to this factor 1/N, the term S 1 is identical to S 1 = ij 2 ij , which is given by the sum of bond length squares 2 ij . The reason why the summation convention (1/N) tet is used instead of the simple form ij for S 1 is that the discrete S 2 is naturally given by the sum over tetrahedrons as in the 2D case described above, and the same summation convention is also assumed for S 1 . The symbol φ i in S 3,i is the internal angle of the triangle (Figure 6c), i in S 3 is the sum over all internal angles i, and σ 1(i) and σ 2(i) in the definition of S 3,i are σ at the bonds {1(i), 2(i)} in which the internal angle is φ i (Figure 6c). The symbol κ for the coefficient of S 3 is the same as the bending rigidity κ for S 3 of the 2D model in Equation (6); however, no confusion is expected. In the case of the 2D model, S 3 only resists pure bending. In contrast, S 3 of the 3D model in Equation (11) resists all deformations of the tetrahedron except the homologous deformation.

Simulation Technique
The standard Metropolis technique is used to update the variables r and σ [21,22]. For the update of r i at vertex i, a new position r i is generated by three uniform random numbers δ r inside a sphere of fixed radius r 0 such that r i = r i +r 0 δ r, where r 0 is a positive number. The new r i is accepted with the probability Min[1, exp −δS], where δS is the change in the Hamiltonian due to the update of r i such that δS = S(new)−S(old) under the constraints U 1 and U 2 on the bond length ij and σ i , respectively. The constraint U 2 is explicitly imposed on the update of σ, and this U 2 also implicitly imposes a constraint on the bond length ij . Indeed, the constraint U 2 on ij in Equation (6) or Equation (11) for σ ij = 1 together with U 1 allows the bond length ij to have only ij = c ij for the rigid bond model and ij < c ij for the empty bond model. The rate of acceptance for the update of r depends on the number r 0 , which is fixed so that the acceptance rate is approximately in the range between 50% and 80%.
The update of σ i to the new σ i is performed independently of the old σ i with the probability Min[1, exp −δS], which is the same expression for the update of r. In this update of σ i , the change from σ i = 0 to σ i = 1, which is called crystallization, is allowed only when i > 0 . This crystallization process, from σ = 0 to σ = 1, causes a sudden decrease in S i (i = 1, 2) for both rigid and empty models. The energy S 3 also discontinuously decreases in the empty model, while it remains unchanged in the rigid bond model. Therefore, the crystallization process is determined only by U i (i = 1, 2) in both rigid and empty bond models. In contrast, the melting process, from σ = 1 to σ = 0, is not constrained by U i (i = 1, 2); therefore, this process is determined only by the energies S i (i = 1, 2, 3) in both rigid and empty bond models. Indeed, in the rigid bond model, a sudden increase is expected in S i (i = 1, 2) for the melting process, and in the case of the empty bond model, a sudden increase is also expected in S i (i = 1, 2, 3) for this process. Therefore, the melting process is determined only by the change in energy in both rigid and empty bond models according to the probability Min[1, exp −δS].
One more point on the update of σ that should be emphasized is that this update is performed only once every 100 Monte Carlo sweeps (MCSs), where one MCS is composed of N consecutive updates of r. The reason why the variables r and σ are updated in such an asymmetrical manner is because a symmetrical or almost symmetrical update of σ causes a configuration with a nonuniform distribution of crystalline bonds, where "symmetrical" means that the total number of updates of r and σ in one MCS is the same. As a consequence of this nonuniform distribution of crystalline bonds, χ for the large strain region becomes lower than that for the intermediate strain region. This effect will be described in further detail in the presentation section.

Frame Tension as Tensile Stress
Detailed information on how to calculate the tensile stress is given in Ref. [20], and, here, we start with the formula. The formula for the stress in the 2D models is given by where k B and T are the Boltzmann constant and the temperature, respectively, A(= πDH) corresponds to the true area of the cylinder, and τ is the simulated frame tension. We should note that 3N − N bnd on the left-hand side should be replaced by 3N − 2N bnd for the 3D model because the integration for the boundary vertices in Z of Equation (2) is the two-dimensional (one-dimensional) integration for the 2D (3D) model. The total number of degrees of freedom 3N is given by 3N = 3N−N cr (rigid bond model) and 3N = 3N−N cr /2 (empty bond model). As mentioned in Section 2.2, the total number of degrees of freedom for the vertex move is slightly reduced by crystallization. For this reason, the number 3N in Equation (13) is replaced by 3N . We briefly describe how to calculate the tensile stress τ sim in Equation (13) and how to change the unit from (N/m) to (Pa). The tensile stress is not directly calculated in our models because the materials represented by 2D or 3D lattices are very thin. Indeed, the 2D lattice is a cylindrical surface, and the 3D lattice is also regarded as a surface because its thickness is almost negligible compared to the diameter. In contrast, the frame tension τ is calculable as a response to a mechanical constraint, which is imposed by fixing the height H of the cylindrical surface (Figure 7a,b). For these reasons, τ is used to obtain τ sim . The calculation formula of τ is actually obtained by a property of the partition function Z in Equation (2) under the scale change such as r → α r (α > 0), which is a change of variables in the multiple-integration of Z. This property in Z is called scale invariance [23], and is expressed by Z(αr) = Z(r) for any α > 0. From this, we have dZ(αr)/dα| α→1 = 0. This derivative of Z is directly calculated from the expression (α r, σ)), which is from the 2D model partition function in Equation (2), where the periodic boundary condition is assumed on the boundary vertices. Hence, the factor of α is calculated by 3(N−N bnd )+2N bnd = 3N−N bnd from the fact that the boundary integration is two-dimensional. Note also that the total number of degrees of freedom 3N should be replaced by 3N , as mentioned above. Except for the total number of degrees of freedom and difference in the symbols for energies, the calculation of dZ(αr)/dα| α→1 is exactly the same as those written in Ref. [20]. Thus, we obtain the surface tension τ = (2 S 1 +4b S 2 −(3N −N bnd ))/2A. By multiplying k B T/a 3 to this τ, we have the stress τ sim in Equation (13) in the unit of (Pa). The reason for the multiplication of k B T/a 3 is as follows: First, we should note that the simulation unit is determined by fixing k B T = 1 (Nm) and a = 1 (m), and, for this reason, all quantities with length units should be multiplied by the lattice spacing a in the physical unit. For this reason, by replacing the area A in the denominator of τ by Aa 2 , we obtain τ/a 2 with the unit of (1/m 2 ). Multiplying k B T to this τ/a 2 , we have τk B T/a 2 with the unit of (N/m). This the physical unit of the frame tension τ. Finally, to obtain the quantity that has the unit of (N/m 2 ), we have to multiply 1/a to this quantity once again, and this leads to the expression of τ sim in Equation (13). By assuming room temperature for T, we have The symbol a(m) is the lattice spacing or the bond length, which corresponds to the coarse-grained distance between polymer segments, and, hence, it can be fixed to an arbitrary number larger than the van der Waals distance (∼10 −10 (m)). This implies that the magnitude of the calculated τ can be controlled arbitrarily such that τ sim is comparable to experimental data τ exp by using this parameter a(> 10 −10 ).
The diameter D, which is equal to the initial height H 0 , is fixed such that the height strain is zero for H = H 0 . For this purpose, MC simulations should be performed to obtain the correct value of H 0 , which depends on the other parameters, such as b and κ, before the start of the production simulations.
Once H 0 is obtained, the next task that should be done before the production simulations is to fix the constant 0 in U 1 of Equation (7). Using the value of H 0 , we fix 0 to 0 = H 0 /10 (16) for both the 2D lattice of size N =10,230 and the 3D lattice of size N = 9760. This value of 0 is crucial to the shape of the crystallization ratio χ vs. strain, which influences the final result of the stress-strain diagram. However, the results are almost independent of a small variation of this value of 0 . The factor 1/10 in 0 = H 0 /10 is expected to be dependent on the size N in general because H 0 depends on L 1 as in Equation (1), and this L 1 determines N such that L 1 = 15 for the 2D lattice of size N = 10,230 and L 1 = 8 for the 3D lattice of size N = 9760. We should note that the value of 0 or H 0 depends on the experimental data of stress-strain curves to be fitted and on whether the model is 2D or 3D. In this paper, we assume Equation (16) for 0 because the experimental data of stress-strain curves are well reproduced, as we see in the following section.

Snapshots
Snapshots for the 2D and 3D empty models are shown in Figure 8a-f, respectively. Small red lines denote the crystalline bonds corresponding to σ = 1. The crystallization ratio χ for the 2D empty model in Figure 8a-c is χ 0.0035, χ 0.060, and χ 0.12, respectively, and χ for the 3D model in Figure 8d-f is χ 0.0016, χ 0.026, and χ 0.083, respectively. The value of χ 0.12 in Figure 8c for the 2D model is close to χ max 0.15, and χ 0.083 in Figure 8f for the 3D model is also close to χ max 0.09. We should note that the crystalline bonds are uniformly distributed in both the outer and inner surfaces in the case of the 3D model.

Results of 2D Models
Pradhan et al. [24] synthesized elastomer nanocomposites with nanoparticle dispersion and studied their mechanical properties. The strain of the carboxylated nitrile rubber (XNBR) with layered double hydroxide (LDH) is reported to be very large (∼1000% or more), and the stress-strain data are influenced by SIC; hence, XNBR/LDH is similar to natural rubber. For this reason, we compare our simulation results with the experimental data of these materials denoted by XL5 and XL10 in Ref. [24]. The numbers 5 and 10 denote the LDH content with the units of phr, parts per hundred rubber.
The reason why these experimental data, reported in [24], are compared with the simulation data is that these reported data are considered to be of an equilibrium state. In contrast, experimental data of natural rubber such as those plotted in Figure 1c are not of equilibrium state, although the data obtained in the unloading process are almost close to those of the equilibrium state.
As described in the introduction, SIC causes a plateau in the intermediate strain region of stress-strain curves, and, moreover, an upturn of the curve or hardening is observed in the large strain region where the plateau terminates. This upturn is also due to the limitation of extensibility of the polymer chains [24]. For this reason, the authors in Ref. [24] used a modified Mooney-Rivlin equation [25,26] to see whether SIC truly influences the stress-strain curves, and found that XL5 and XL10 are influenced by SIC. Here, we should note that, as mentioned in the introduction, we focus on the stress relaxation in this paper. For this reason, the upturn in the large strain region is simply implemented in the models by the term S 2 , which is quadratic with respect to the bond length ij , and, therefore, it is almost independent of crystalline bonds.
In Figure 9a, the experimental data (×:Exp) of XL5 and the simulation results ( and ) of the 2D rigid bond model are plotted. The parameters b and κ are fixed to b = 0.003 and κ = 0.85, respectively. The coefficient b of S 2 , which is the quadratic part of tensile energy, is fixed to a nonzero value, though it is very small compared to the coefficient 1 of S 1 . This small but nonzero b plays a role in the large strain region to increase τ sim in both models, as mentioned above. The assumed value of 0 in Equation (16) for 2D models is given by The simulation data ( ) correspond to χ = 0, which implies that no crystalline bond is included. We find that no plateau region is included in the data of χ = 0. In contrast, the curve of data ( ) is clearly different from ( ) and has a small plateau at the region 2 < ε < 4, where χ plotted in Figure 9b increases from χ = 0 to χ max ( 0.15). This result implies that τ decreases with increasing χ, and, hence, this corresponds to the stress relaxation induced by SIC. The reason why τ ( ) becomes lower than τ ( ) is very simple to understand because the crystalline bond is defined to have no tensile energy, and, therefore, S 1 , S 2 , and 3N in τ sim of Equation (13) decrease with increasing χ. The results of the 2D empty bond model are plotted in Figure 9c,d, which are almost the same as those of the 2D rigid bond model shown in Figure 9a,b. The parameters b and κ in the 2D empty model are fixed to b = 0.002 and κ = 0.85, respectively, which are almost the same as those fixed in the 2D rigid model.
Thus, the stress relaxation is understood to be a result of the fact that the total tensile energy accumulated in the amorphous part is simply decreased by crystallization in real materials. Indeed, the total number of segments of the amorphous part that shares the tensile energy is decreased by SIC in the models. Moreover, in real materials, at the small strain region where the crystallization starts to increase, almost no accumulation of the tensile energy is expected in the crystalline part. This property in real materials is shared by the models in this paper. However, the tensile energy shared by the crystalline part cannot be neglected for the large strain region, where the amorphous part is almost maximally extended and has no room for the accumulation of tensile energy. This property is not shared by the models in this paper, and, for this reason, the sudden upturn by SIC at the large strain region is not reproduced. We consider that a part of the difference between the curves ( ) and ( ) corresponds to the entropy change. The temperature in the recovery process is also expected to be lower than that in the extension process in real experiments, and this temperature difference is due to the decrease in entropy from the thermodynamics viewpoint. The decrease in entropy from the microscopic viewpoint is naturally expected from our models because the crystalline bond of σ = 1 effectively decreases the total number of degrees of freedom for the variables r i , as mentioned in Section 2.2.
In Figure 10a-f, χ vs. MCS are plotted, and the convergence of MC simulations for 2D rigid and empty bond models are shown. We find that in both models, the convergence at the transition region (Figure 10b,e), where χ has an intermediate value between 0 and χ max , is relatively slower than those at small (Figure 10a,d) and large strain regions (Figure 10c,f). The slow speed of convergence at the intermediate strain region is reasonable because crystallization is understood to be a first-order transition. One more point to note is that χ has an intermediate value between χ = 0 and χ max in Figure 10b,e. This observation indicates that χ smoothly changes with respect to strain ε even though SIC is a first-order transition. We will not go into detail on this first-order transition because our focus in this paper is not on the first-order nature of the crystallization.
Finally, in this subsection, we should comment on the reason for why the MC update of σ is performed only once every 100 MCSs. As mentioned in Section 2.4, the crystallization ratio χ at the large strain region becomes slightly smaller than χ max if the MC update is performed more frequently. This result comes from the fact that the convergence speed of the variable r is very slow compared with that of σ. In other words, if σ is updated more frequently before the lattice configuration is equilibrated, then the crystalline bonds of σ = 1 are expected to be distributed non-uniformly, and χ will be smaller than χ max at the large strain region.

Results of the 3D Model
Now, we show the simulation results of the 3D empty bond model. The main difference between 3D and 2D models is that the Hamiltonian depends on the dimension D, and this difference in dimension is a huge difference in general from the modeling view point, and, therefore, it is interesting to see whether or how the results change depending on the dimension D. The reason for why only the empty bond model is simulated in the 3D case is because no specific difference is found between the two models, rigid and empty, in the 2D case. The obtained results τ sim and χ are plotted in Figure 11a-d together with the reported experimental data. The symbols (×:Exp) in Figure 11a,c correspond to the experimental data denoted by XL5 and XL10, respectively [24]. The parameters (b, κ) are fixed to (b, κ) = (0.0015, 0.3) and (b, κ) = (0.003, 0) for the experimental data XL5 and XL10, respectively. The value of 0 in Equation (16) for the 3D model is given by 0 = 2.69a (3D empty, XL5), 0 = 2.84a (3D empty, XL10). (18) We find that the simulation results by the 3D empty model are in good agreement with the experimental data, just as in the case of the 2D models.
We plot χ vs. MCS obtained at the small, intermediate, and large strain regions in Figure 12a-c to see the convergence of the MC simulations. It is clear that the convergence speed in the intermediate strain region is also slower than that in the small and large strain regions, as in the 2D models. Note also that the convergence speed of the 3D model is relatively faster than that of the 2D models. We consider that this result comes from the fact that the surface or lattice fluctuation of the 3D models is relatively smaller than that of the 2D model, or, in other words, the phase space volume for the MC update of the polymer position r i is expected to be relatively smaller in the 3D model than in the 2D models. Indeed, the 2D lattice is composed of triangles, and, hence, the vertex position easily moves even when κ is fixed to κ = 0.85. In contrast, the 3D lattice is composed of tetrahedrons, and, hence, the vertices hardly move compared to those on the 2D lattice.  Next, we calculate the assumed lattice spacing a to calculate τ, as plotted in Figures 9 and 10. The results are shown in Table 1. The lattice spacing is considered to be the unit of distance between two neighboring vertices, which are understood to be coarse-grained lengths of polymer segments (see Equations (17) and (18)). Therefore, the value of a should be at least larger than the van der Waals distance ( 10 −10 [m]). From Table 1, we find that all assumed a values are larger than 10 −10 [m], and, therefore, lattice modeling is meaningful as a coarse-grained technique for polymer networks. Table 1. The lattice spacing a used for τ sim (in Equation (14)) plotted in Figures 9 and 11. Finally, in this subsection, we plot the mean length cr for crystalline bonds in Figure 13a,b. This cr is not always identical to c i in Equation (7) because of the constraint cr < c i for the crystalline bond length in the empty bond model. However, as described in Section 2.2, c i is dynamically changeable, and, hence, cr is also dynamically changeable and depends on the strain ε. The sudden jump of c i in Figure 13a for the 2D models implies that there are no crystalline bonds in the small strain region close to ε 0, and the nonzero c i in Figure 13b implies that there is more than one crystalline bond even in the small strain region in the case of the 3D model. We find that c i increases with increasing ε, and this behavior is consistent with that observed in real materials [3].

Dependence on Stiffness κ
In the preceding subsection, two different experimental data, XL5 and XL10, are simulated and plotted in Figure 11a,c, where the stiffness κ = 0.3 and κ = 0 are assumed, respectively. However, the dependence of the results on κ is not so clear. Therefore, in this subsection, we should discuss the effect of stiffness κ on the stress-strain curve and clarify that the effect of κ is different from the effect of SIC. First, we plot in Figure 14a the results of 3D model with the crystalline state ( ) (⇔ χ > 0) and those without the crystalline state ( , ) (⇔ χ = 0), where κ is fixed to κ = 1 and κ = 0, respectively. In Figure 14a, the values of τ are given by Equation (13) without the factor k B T/a 3 , and these are raw simulation data expressed by the simulation unit. It is clear from Figure 14a that the data ( ) decrease in the intermediate strain region and are smaller than the other two data ( , ) of χ = 0, except in the smaller strain region. This decrement of τ ( ) is understood to be due to the crystalline states.
To see the effect of κ more clearly, the stress τ sim , which includes the factor k B T/a 3 in Equation (13), is plotted in Figure 14b with the experimental data (×). Inside Figure 14b, the data in the small strain region are drawn by solid and dashed lines. We find from these solid and dashed lines that the increment of κ increases the stress only at the region ε ≤ 2.5. On the other hand, SIC has no effect on the stress in the small strain region in general, because SIC is effective only at the large strain region, such as ε ≥ 2.5, for example. Thus, we find from these numerical data that the stress relaxation by SIC is independent of the effect of κ. In other words, the stress relaxation by SIC, observed in experimental data, is unable to be understood without the crystalline states, such as the rigid or empty bond introduced in Section 2.2, at least in the framework of standard Gaussian chain models.
We also find that the results of the 2D models without crystalline states are completely inconsistent with the experimental data even when κ is increased to be sufficiently large. This means that, in the case of 2D models, experimental data can only be reproduced with the crystalline states, which are introduced by rigid bonds or empty bonds in Equations (4) Figure 14. (a) The stress τ vs. ε obtained by the 3D model with the crystalline state ( ) and without the crystalline state ( , ), and (b) the stress τ sim vs. ε obtained without the crystalline state. The stress τ in (a) has the simulation unit, which corresponds to k B T = 1 and is not modified by the factor k B T/a 3 in Equation (13), while the stress τ sim in (b) has the real physical unit (MPa) and is modified by the factor k B T/a 3 . Experimental data denoted by Exp (×) in (b) are XL5, the same as in Figure 11c.

Summary and Conclusions
We have studied the equilibrium property of stress relaxation observed in rubbers undergoing strain-induced crystallization (SIC) by Monte Carlo (MC) simulations for simple models constructed by extending the Gaussian chain model. Two types of new models, rigid and empty, are introduced and studied. In these models, a new variable σ is introduced, and it is defined on the bond to distinguish the amorphous and crystalline states. The crystalline state is simply defined to have no tensile energy, where the tensile energy is given by a spring potential or the so-called Gaussian bond potential and is defined only on the amorphous state.
We find that the obtained MC results for the stress-strain curves are in good agreement with the existing experimental data of composite materials, of which the mechanical property was reported to be close to that of natural rubber, and the corresponding stress-strain curves are influenced by SIC [24]. Here, we have to emphasize that these experimental data are not time-dependent ones, which can be considered the equilibrium ones, and, hence, can be close to unloading stress-strain curves. In the simulation data, stress relaxation is clearly observed in the intermediate strain region, where the crystallization ratio χ starts to increase. The fact that no tensile energy is accumulated in the crystalline bonds in our models indicates that the stress relaxation is simply due to the decrease in the tensile energy of the amorphous state, where this decrease in energy is in the sense that the total energy is decreased due to the change in state from amorphous to crystalline. No contribution to the tensile energy is expected from the crystalline state, at least in the intermediate strain region.
In this paper, we focus only on stress relaxation, and the upturn of the stress-strain curve or hardening at the large strain region is implemented simply by including a quadratic term in the tensile energy, which is not directly connected to the crystalline state. This upturn of the stress-strain curve is another interesting phenomenon associated with SIC, and, hence, it remains to be studied more intuitively from a microscopic viewpoint.

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

Abbreviations
The following abbreviations are used in this manuscript:

2D
two-dimensional 3D three-dimensional MC Monte Carlo MCS Monte Carlo Sweep SIC Strain-Induced Crystallization