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

Vladislav Egorov, Hiroshi Koibuchi, Chrystelle Bernard, Jean-Marc Chenal, Gildas Diguet, Gaël Sebald, Jean-Yves Cavaille, Toshiyuki Takagi, Laurent Chazeau

► To cite this version:

Vladislav Egorov, Hiroshi Koibuchi, Chrystelle Bernard, Jean-Marc Chenal, Gildas Diguet, et al. Coarse-Grained Lattice Modeling and Monte Carlo Simulations of Stress Relaxation in Strain-Induced Crystallization of Rubbers. Polymers, 2020, 12 (6), pp.1267. 10.3390/polym12061267. hal-02970928

HAL Id: hal-02970928

https://hal.science/hal-02970928v1

Submitted on 14 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

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

Vladislav Egorov <fnref n="1" />,*, Hiroshi Koibuchi <fnref n="2" />,*, Chrystelle Bernard <fnref n="3" />,<fnref n="4" />, Jean-Marc Chenal <fnref n="5" />,
Gildas Diguet <fnref n="3" />,<fnref n="5" />, Gael Sebald <fnref n="3" />,<fnref n="5" />, Jean-Yves Cavaille <fnref n="3" />,<fnref n="5" />, Toshiyuki Takagi <fnref n="6" />
and Laurent Chazeau <fnref n="5" />

<figure><img src="image_2.png" /></figure>

check for
updates

Received: 7 May 2020; Accepted: 27 May 2020; Published: 1 June 2020

**Abstract:** 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.

**Keywords:** rubber elasticity; strain-induced crystallization; stress relaxation; Monte Carlo; stress–strain curves; statistical mechanics

## 1. Introduction

Natural rubbers undergo strain-induced crystallization (SIC), which attracts a lot of attention and has been studied extensively [1–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

<footnote n="1">¹ Cherepovets State University (ChSU), Prospekt Lunacharskogo, 5, 162600 Cherepovets, Vologda Oblast, Russia; rvladegorov@rambler.ru</footnote>

<footnote n="2">² National Institute of Technology (KOSEN), Sendai College, 48 Nodayama, Medeshima-Shiote, Natori-shi, Miyagi 981-1239, Japan</footnote>

<footnote n="3">³ ELyTMAx UMI 3757, CNRS-Université de Lyon, Tohoku University, International Joint Unit, Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai 980-8577, Japan; chrystelle.bernard@rift.mech.tohoku.ac.jp (C.B.); gildas.diguet.d4@tohoku.ac.jp (G.D.); gael.sebald@insa-lyon.fr (G.S.); jean-yves.cavaille@insa-lyon.fr (J.-Y.C.)</footnote>

<footnote n="4">⁴ Frontier Research Institute for Interdisciplinary Sciences (FRIS), Tohoku University, 6-3 Aoba Aramaki, Aoba-ku, Sendai 980-8578, Japan</footnote>

<footnote n="5">⁵ Materials Engineering and Science (MATEIS), CNRS, INSA Lyon UMR 5510, Université de Lyon Batiment B. Pascal, Avenue Jean Capelle, CEDEX, 69621 Villeurbanne, France; jean-marc.chenal@insa-lyon.fr (J.-M.C.); laurent.chazeau@insa-lyon.fr (L.C.)</footnote>

<footnote n="6">⁶ Tohoku Forum for Creativity, Tohoku University, 2-1-1 Katahira, Aoba-ku Sendai 980-8577, Japan; toshiyuki.takagi.d4@tohoku.ac.jp</footnote>

<footnote>* Correspondence: koi-hiro@sendai-nct.ac.jp or koibuchih@gmail.com</footnote>
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.

<figure><img src="image_3.png" /><figcaption>**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 $\chi$ 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 $\tau$ (MPa) and the crystallization ratio $\chi$ (%) in (c) are taken from experimental data in [5].</figcaption></figure>

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–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–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.

## 2. Models

### 2.1. 2D Lattice

A 2D triangulated lattice was originally constructed for the surface model of membranes [14–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

$$H_0 = D = (\sqrt{3}/2)(L_1 - 1)a, \quad (1)$$

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 $\pi D = (L_2-1)a$. Therefore, from $D = H_0$, the relation $L_2 = (\sqrt{3}\pi/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 = (\sqrt{3}\pi/2)(L_1-1)+1$ is not always exactly satisfied because the right-hand side is not an integer.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) A cylindrical surface is made of a rectangular surface of size $L_1 \times 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$.</figcaption></figure>
## 2.2. Rigid Bond Model and Empty Bond Model

Let $\vec{r}_i \in \mathbf{R}^3$, ($i=1, \dots, N$) be the vertex position of the lattice (Figure 3a), $\ell_{ij} = |\vec{r}_i - \vec{r}_j|$ be the bond length (Figure 3b), and $\vec{n}_i$ be the unit normal vector of the triangle (Figure 3c). These variables $\vec{r}_i$ are integrated into the partition function $Z$ for 2D models (here, we discuss 2D surface models), and $Z$ is given by

$$Z(\ell_0) = \sum_{\sigma} \int \prod_{i=1}^{N-N_{\text{bnd}}} d\vec{r}_i \prod_{i=1}^{N_{\text{bnd}}} d\vec{r}_i \exp(-S(\vec{r}, \sigma)), \quad (2)$$

where the symbol $\ell_0$ in $Z(\ell_0)$ indicates that the model implicitly depends on a constant $\ell_0$, which will be used to define the crystalline bond below. The symbol $\int \prod_{i=1}^{N-N_{\text{bnd}}} d\vec{r}_i$ denotes the 3D integration for the internal vertices. The other $\int \prod_{i=1}^{N_{\text{bnd}}} d\vec{r}_i$ denotes 2D integration for the boundary vertices on the upper and lower boundaries with a periodic boundary condition such that

$$Z_i^U = Z_i^L + H, \quad X_i^2 + Y_i^2 = D^2/4, \quad (3)$$

where $Z_i^U$ and $Z_i^L$ denote the Z component of the vertex position $\vec{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 $\int \prod_{i=1}^{N_{\text{bnd}}} d\vec{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, $\int \prod_{i=1}^{N_{\text{bnd}}} d\vec{r}_i$ becomes a one-dimensional integration because no periodic boundary condition is imposed.

<figure><img src="image_5.png" /><figcaption>Figure 3. (a) A triangular mesh or lattice and the vertex position $\vec{r}_i$, (b) a triangle with vertices 1, 2, 3 and edge lengths $\ell_{12}$, $\ell_{23}$, $\ell_{31}$, and (c) the unit normal vectors $\vec{n}_0$ of the triangle 123, and those $\vec{n}_i$, ($i=1, 2, 3$) of its neighboring triangles, where the suffix i denotes the corresponding triangle.</figcaption></figure>

The symbol $\sigma$ in Equation (2), which has values in $\{0, 1\}$, is another dynamical variable defined on the bond $i$. The value of $\sigma$ represents whether the bond $i$ is rigid or amorphous:

$$\sigma_i = \begin{cases} 0 & (\Leftrightarrow \text{bond } i = \text{amorphous}) \\ 1 & (\Leftrightarrow \text{bond } i = \text{rigid}) \end{cases}, \quad (\text{rigid bond model}). \quad (4)$$

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 $\ell_{ij} (=|\vec{r}_i - \vec{r}_j|) = \ell_{ij}^c$, where $\ell_{ij}^c$ is described in the following paragraph. In contrast, this constraint on $\ell_{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

$$\sigma_i = \begin{cases} 0 & (\Leftrightarrow \text{bond } i = \text{amorphous}) \\ 1 & (\Leftrightarrow \text{bond } i = \text{empty}) \end{cases}, \quad (\text{empty bond model}). \quad (5)$$
The variable $\sigma_i$ in Equations (4) and (5) can be changed from $\sigma_i=0$ to $\sigma_i=1$ when the bond length $\ell_i$ is $\ell_i > \ell_0$ using the constant $\ell_0$ in $Z(\ell_0)$ according to the Monte Carlo update procedure described in a later section. If $\sigma_i$ is changed from $\sigma_i=0$ to $\sigma_i=1$, then the bond length $\ell_i$ is fixed to the length $\ell_i^c$, which is the length of bond $i$ just before the change in $\sigma_i$ in the rigid model. In the case of the empty bond model, the length $\ell_i$ of the crystalline bond is changeable to a new length $\ell'_i$ if $\ell'_i < \ell_i^c$, and this $\ell_i^c$ is also changeable depending on the update of $\sigma_i$, which is the same as in the case of the rigid bond model. The constraint $\ell_i^c > \ell_0$ is always satisfied in both models. The empty bond is almost close to the rigid bond except that the length $\ell_i$ is changeable to $\ell'_i < \ell_i^c$. This fact implies that the empty bond length $\ell_i$ is also prohibited from extending to $\ell'_i > \ell_i^c$. The bond configurations for amorphous, rigid, and empty bonds are illustrated in Figure 4a-c.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Illustration of the bond $i$ for (a) $\sigma_i = 0$ corresponding to the amorphous bond, (b) $\sigma_i = 1$ corresponding to the rigid bond, and (c) $\sigma_i = 1$ corresponding to the empty bond.</figcaption></figure>

We should comment that the total number of degrees of freedom of variable $\vec{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 $\vec{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

$$ 
\begin{aligned}
S &= S_1 + bS_2 + \kappa S_3 + U_1 + U_2, \\
S_1 &= \frac{1}{2} \sum_{\Delta} S_{1\Delta}, \quad S_{1\Delta} = \ell_{12}^2 + \ell_{23}^2 + \ell_{31}^2, \\
S_2 &= \frac{1}{2} \sum_{\Delta} (S_{2\Delta}^{(1)} + S_{2\Delta}^{(2)}), \quad S_{2\Delta}^{(1)} = \ell_{12}^4 + \ell_{23}^4 + \ell_{31}^4, \quad S_{2\Delta}^{(2)} = 2(\ell_{12}^2\ell_{13}^2 + \ell_{21}^2\ell_{23}^2 + \ell_{32}^2\ell_{31}^2), \\
\ell_{ij} &= \begin{cases} |\vec{r}_j - \vec{r}_i| & (\sigma = 0) \\ 0 & (\sigma = 1) \end{cases}, \\
S_3 &= \frac{1}{2} \sum_{\Delta} (S_{3\Delta,1} + S_{3\Delta,2} + S_{3\Delta,3}), \quad S_{3\Delta,i} = \begin{cases} 1 - \vec{n}_0 \cdot \vec{n}_i & (\sigma = 0) \\ 1 - \vec{n}_0 \cdot \vec{n}_i \text{ (rigid) or 0 (empty)} & (\sigma = 1) \end{cases}
\end{aligned} 
\quad (6) 
$$
and

$$
\begin{equation}
\begin{aligned}
U_1 &= \sum_{\text{bond } i} U_{1i}(\vec{r}, \sigma), \quad U_{1i}(\vec{r}, \sigma) = \begin{cases} \infty & (\sigma_i = 1 \ \& \ \ell_i \neq \ell_i^c \text{ : rigid}) \text{ or } (\sigma_i = 1 \ \& \ \ell_i > \ell_i^c \text{ : empty}) \\ 0 & (\text{otherwise}) \end{cases}, \\
U_2 &= \sum_{\text{bond } i} U_{2i}(\vec{r}, \sigma), \quad U_{2i}(\vec{r}, \sigma) = \begin{cases} \infty & (\sigma_i = \sigma_j = 1 \ \& \ i \text{ and } j \text{ are connected}) \\ 0 & (\text{otherwise}) \end{cases},
\end{aligned}
\end{equation}
$$

where the first two terms $S_1+bS_2$ and the third term $\kappa 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 $\sigma$ in the definition of $\ell_{ij}$ is the variable $\sigma$ on the bond connecting the vertices $i$ and $j$. No difference between the rigid and empty bond models is in the definition of $\ell_{ij}$ except that the length $\ell_{ij}$ of crystalline bond $ij$ is allowed to be $\ell_{ij} < \ell_{ij}^c$ in the case of the empty bond model. In contrast, the length of bond $ij$ is prohibited from being $\ell_{ij} > \ell_{ij}^c$ for $\sigma_{ij}=1$ in both models. In this sense, this definition of $\ell_{ij}$ is specific to the models for SIC in this paper. In the bending energy $S_{3\Delta,i}(i=1,2,3)$, $\sigma$ 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 $\ell_{ij}$ and $1-\vec{n}_i \cdot \vec{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 = \sum_{\text{bond } i} \ell_i^2$, which is given by the sum of bond length squares $\ell_i^2$. The reason why the sum $\sum_{\Delta}$ is used for $S_1$ instead of $\sum_{\text{bond } i}$ is because the discrete form of the quadratic term $S_2$ is naturally given by $\sum_{\Delta}$, 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 $\ell_i$ in both models. From this constraint $U_1$, the length $\ell_i$ of the empty bond is only allowed to be $\ell_i \le \ell_i^c$, as mentioned above; however, this is strictly prohibited for the rigid bond because $\ell_i$ should always satisfy $\ell_i = \ell_i^c$ from the definition of $\ell_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 $\chi$ such that $\chi \le \chi_{\text{max}}$, where $\chi$ is defined by

$$ \chi = (1/N_B) \sum_{\text{bond } i} \sigma_i, \qquad (8) $$

where $\sum_i \sigma_i$ is the total number of crystalline bonds and $N_B(=\sum_i 1)$ is the total number of bonds. The value of $\chi_{\text{max}}$ depends on the dimension of the lattice, that is, whether it is a 2D or 3D lattice; however, in both cases, $\chi_{\text{max}}$ is reasonable in the sense that it is close to the experimental value, such as $\chi_{\text{max}}(\text{Exp}) \approx 0.1$ [1–4], which will be shown in the presentation section.

The mean value $\langle Q \rangle$ of physical quantity $Q(\vec{r}, \sigma)$ (where the symbol $Q$ denotes an arbitrary quantity) is given by

$$ \langle Q \rangle = \sum_{\sigma} \int \prod_{i=1}^{N-N_{\text{bnd}}} d\vec{r}_i \prod_{i=1}^{N_{\text{bnd}}} d\vec{r}_i Q(\vec{r}, \sigma) \exp(-S(\vec{r}, \sigma))/Z. \qquad (9) $$

We simply write the mean value $Q$ by removing the symbol $\langle \cdot \rangle$ henceforth. To calculate these multiple integrations and the sum $\sum_{\sigma}$ over all possible configurations of $\sigma$, 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 $\vec{r}$ and $\sigma$, which are the vertex position $\vec{r} = \{\vec{r}_1, \cdots, \vec{r}_N\}$ and the spin variable at the vertices $\sigma = \{\sigma_1, \cdots, \sigma_N\}$, respectively.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) A configuration that is not allowed because the two bonds of $\sigma=1$ are connected, and (b) a configuration that is allowed because the two bonds of $\sigma=1$ are separated. The maximal value $\chi_{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.</figcaption></figure>

The 2D continuous forms of the Hamiltonian corresponding to the discrete $S_1$ and $S_2$ in Equation (6) are given by

$$S_1 = \int \sqrt{g} d^2 x g^{ab} \frac{\partial \vec{r}}{\partial x^a} \cdot \frac{\partial \vec{r}}{\partial x^b}, \quad S_2 = \int \sqrt{g} d^2 x \left( g^{ab} \frac{\partial \vec{r}}{\partial x^a} \cdot \frac{\partial \vec{r}}{\partial x^b} \right)^2, \qquad (10)$$

where $g^{ab}$ is the inverse of a metric tensor $g_{ab}$ given by the $2 \times 2$ matrix, and $g$ is the determinant. To obtain the discrete Hamiltonian, we simply assume that $g_{ab} = \delta_{ab}$, which is called the Euclidean metric. By the replacements $\partial_1 \vec{r} (= \partial \vec{r} / \partial x^1) \to \vec{r}_2 - \vec{r}_1$ and $\partial_2 \vec{r} \to \vec{r}_3 - \vec{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].

## 2.3. 3D Model

The three-dimensional (3D) continuous Hamiltonian is obtained from the 2D version in Equation (10) simply by replacing the 2D integration $\int \sqrt{g} d^2 x$ with the 3D integration $\int \sqrt{g} d^3 x$. In the 3D case, the metric tensor $g_{ab}$ is also replaced by a $3 \times 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.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) A cylinder for the 3D model, (b) a tetrahedron with vertices 1, 2, 3, 4 and the corresponding edge length $\ell_{ij}$, ($i \neq j = 1, \dots, 4$), and (c) an internal angle $\phi_1$ of the triangle 123, where 1(1) and 2(1) are the edges, of which the internal angle is $\phi_1$.</figcaption></figure>
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.

$$ 
\begin{aligned}
S &= S_1 + bS_2 + \kappa S_3 + U_1 + U_2, \\
S_1 &= \frac{1}{\bar{N}} \sum_{\text{tet}} S_{1\text{tet}}, \quad S_{1\text{tet}} = \ell_{12}^2 + \ell_{13}^2 + \ell_{14}^2 + \ell_{23}^2 + \ell_{24}^2 + \ell_{34}^2, \\
S_2 &= \frac{1}{\bar{N}} \sum_{\text{tet}} (S_{2\text{tet}}^{(1)} + S_{2\text{tet}}^{(2)}), \quad S_{2\text{tet}}^{(1)} = \ell_{12}^4 + \ell_{13}^4 + \ell_{14}^4 + \ell_{23}^4 + \ell_{24}^4 + \ell_{34}^4, \\
S_{2\text{tet}}^{(2)} &= 2(\ell_{12}^2\ell_{13}^2 + \ell_{12}^2\ell_{14}^2 + \ell_{13}^2\ell_{14}^2 + \ell_{21}^2\ell_{23}^2 + \ell_{21}^2\ell_{24}^2 + \ell_{23}^2\ell_{24}^2 \\
&\quad + \ell_{31}^2\ell_{32}^2 + \ell_{31}^2\ell_{34}^2 + \ell_{32}^2\ell_{34}^2 + \ell_{41}^2\ell_{42}^2 + \ell_{41}^2\ell_{43}^2 + \ell_{42}^2\ell_{43}^2),
\end{aligned} 
\quad (11) $$

$$ \ell_{ij} = \begin{cases} |\vec{r}_j - \vec{r}_i| & (\sigma = 0) \\ 0 & (\sigma = 1) \end{cases}, $$

$$ S_3 = \sum_{\text{triangle } i} S_{3,i}, $$

$$ S_{3,i} = \begin{cases} 1 - \cos(\phi_i - \pi/3) & (\sigma_{1(i)} = 0 \text{ and } \sigma_{2(i)} = 0, \text{ or empty}) \\ 0 & (\text{otherwise}) \end{cases}, $$

where $\sum_{\text{tet}}$ in $S_1$ and $S_2$ denotes the sum over the tetrahedrons (Figure 6a). The symbol $\bar{N}$ is the mean value of the total number of tetrahedrons around a bond and is given by

$$ \bar{N} = (1/N_B) \sum_{\text{bond } i} \sum_{\text{tet}(i)} 1, \quad (12) $$

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/\bar{N}$, the term $S_1$ is identical to $S_1 = \sum_{ij} \ell_{ij}^2$, which is given by the sum of bond length squares $\ell_{ij}^2$. The reason why the summation convention $(1/\bar{N}) \sum_{\text{tet}}$ is used instead of the simple form $\sum_{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 $\phi_i$ in $S_{3,i}$ is the internal angle of the triangle (Figure 6c), $\sum_i$ in $S_3$ is the sum over all internal angles $i$, and $\sigma_{1(i)}$ and $\sigma_{2(i)}$ in the definition of $S_{3,i}$ are $\sigma$ at the bonds $\{1(i), 2(i)\}$ in which the internal angle is $\phi_i$ (Figure 6c). The symbol $\kappa$ for the coefficient of $S_3$ is the same as the bending rigidity $\kappa$ 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.

## 2.4. Simulation Technique

The standard Metropolis technique is used to update the variables $\vec{r}$ and $\sigma$ [21,22]. For the update of $\vec{r}_i$ at vertex $i$, a new position $\vec{r}_i'$ is generated by three uniform random numbers $\delta\vec{r}$ inside a sphere of fixed radius $r_0$ such that $\vec{r}_i' = \vec{r}_i + r_0\delta\vec{r}$, where $r_0$ is a positive number. The new $\vec{r}_i'$ is accepted with the probability Min[1, exp $-$ $\delta S$], where $\delta S$ is the change in the Hamiltonian due to the update of $\vec{r}_i$ such that $\delta S = S(\text{new}) - S(\text{old})$ under the constraints $U_1$ and $U_2$ on the bond length $\ell_{ij}$ and $\sigma_i$, respectively. The constraint $U_2$ is explicitly imposed on the update of $\sigma$, and this $U_2$ also implicitly imposes a constraint on the bond length $\ell_{ij}$. Indeed, the constraint $U_2$ on $\ell_{ij}$ in Equation (6) or Equation (11) for $\sigma_{ij} = 1$ together with $U_1$ allows the bond length $\ell_{ij}$ to have only $\ell_{ij} = \ell_{ij}^c$ for the rigid bond model and $\ell_{ij} < \ell_{ij}^c$ for the empty bond model. The rate of acceptance for the update of $\vec{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 $\sigma_i$ to the new $\sigma'_i$ is performed independently of the old $\sigma_i$ with the probability Min[1, exp $-$ $\delta S$], which is the same expression for the update of $\vec{r}$. In this update of $\sigma_i$, the change
from $\sigma_i=0$ to $\sigma_i=1$, which is called crystallization, is allowed only when $\ell_i > \ell_0$. This crystallization process, from $\sigma=0$ to $\sigma=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 $\sigma=1$ to $\sigma=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 $-\delta S$].

One more point on the update of $\sigma$ 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 $\vec{r}$. The reason why the variables $\vec{r}$ and $\sigma$ are updated in such an asymmetrical manner is because a symmetrical or almost symmetrical update of $\sigma$ causes a configuration with a nonuniform distribution of crystalline bonds, where “symmetrical” means that the total number of updates of $\vec{r}$ and $\sigma$ in one MCS is the same. As a consequence of this nonuniform distribution of crystalline bonds, $\chi$ 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.

## 2.5. 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

$$ \tau_{\text{sim}} = \frac{2\langle S_1 \rangle + 4b\langle S_2 \rangle - (3N' - N_{\text{bnd}}) \frac{k_B T}{a^3}}{2A} = \tau \frac{k_B T}{a^3} \quad (\text{Pa}), \quad (13) $$

where $k_B$ and $T$ are the Boltzmann constant and the temperature, respectively, $A(=\pi DH)$ corresponds to the true area of the cylinder, and $\tau$ is the simulated frame tension. We should note that $3N' - N_{\text{bnd}}$ on the left-hand side should be replaced by $3N' - 2N_{\text{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_{\text{cr}}$ (rigid bond model) and $3N' = 3N - N_{\text{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 $\tau_{\text{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 $\tau$ 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, $\tau$ is used to obtain $\tau_{\text{sim}}$. The calculation formula of $\tau$ is actually obtained by a property of the partition function $Z$ in Equation (2) under the scale change such as $\vec{r} \to \alpha\vec{r}$ ($\alpha > 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(\alpha\mathbf{r}) = Z(\mathbf{r})$ for any $\alpha > 0$. From this, we have $dZ(\alpha\mathbf{r})/d\alpha|_{\alpha \to 1} = 0$. This derivative of $Z$ is directly calculated from the expression $Z(\alpha\mathbf{r}) = \alpha^{3N-N_{\text{bnd}}} \sum_{\sigma} \int \prod_{i=1}^{N-N_{\text{bnd}}} d\vec{r}_i \prod_{i=1}^{N_{\text{bnd}}} d\vec{r}_i \exp(-S(\alpha\vec{r}, \sigma))$, 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 $\alpha$ is calculated by $3(N-N_{\text{bnd}})+2N_{\text{bnd}}=3N-N_{\text{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(\alpha\mathbf{r})/d\alpha|_{\alpha \to 1}$ is exactly the same as those written in Ref. [20]. Thus, we obtain the surface tension $\tau = (2\langle S_1 \rangle + 4b\langle S_2 \rangle - (3N' - N_{\text{bnd}}))/2A$. By multiplying
$k_B T / a^3$ to this $\tau$, we have the stress $\tau_{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 $\tau$ by $Aa^2$, we obtain $\tau / a^2$ with the unit of $(1/m^2)$. Multiplying $k_B T$ to this $\tau / a^2$, we have $\tau k_B T / a^2$ with the unit of (N/m). This the physical unit of the frame tension $\tau$. 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 $\tau_{sim}$ in Equation (13).

<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Frame or projected area to be fixed to measure the surface tension $\tau$, which can be called frame tension, in the case of a square surface, and (b) projected area $A = \pi DH$ is fixed to calculate $\tau$ in the case of a cylindrical surface.</figcaption></figure>

By assuming room temperature for $T$, we have

$$ \tau_{sim} = (4 \times 10^{-21}) \frac{\tau}{a^3} \text{ (Pa)}. \qquad (14) $$

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 $\tau$ can be controlled arbitrarily such that $\tau_{sim}$ is comparable to experimental data $\tau_{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

$$ \varepsilon = H/H_0 - 1 \qquad (15) $$

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 $\kappa$, 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 $\ell_0$ in $U_1$ of Equation (7). Using the value of $H_0$, we fix $\ell_0$ to

$$ \ell_0 = H_0 / 10 \qquad (16) $$

for both the 2D lattice of size $N = 10,230$ and the 3D lattice of size $N = 9760$. This value of $\ell_0$ is crucial to the shape of the crystallization ratio $\chi$ 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 $\ell_0$. The factor $1/10$ in $\ell_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 $\ell_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 $\ell_0$ because the experimental data of stress-strain curves are well reproduced, as we see in the following section.
## 3. Monte Carlo Results

### 3.1. 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 $\sigma=1$. The crystallization ratio $\chi$ for the 2D empty model in Figure 8a–c is $\chi \approx 0.0035$, $\chi \approx 0.060$, and $\chi \approx 0.12$, respectively, and $\chi$ for the 3D model in Figure 8d–f is $\chi \approx 0.0016$, $\chi \approx 0.026$, and $\chi \approx 0.083$, respectively. The value of $\chi \approx 0.12$ in Figure 8c for the 2D model is close to $\chi_{max} \approx 0.15$, and $\chi \approx 0.083$ in Figure 8f for the 3D model is also close to $\chi_{max} \approx 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.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Snapshots of the (a–c) 2D empty model and (d–f) 3D model, where small red lines denote the crystalline bonds. The crystallization ratio is (a) $\chi \approx 0.0035$, (b) $\chi \approx 0.060$, (c) $\chi \approx 0.12$, (d) $\chi \approx 0.0016$, (e) $\chi \approx 0.026$, and (f) $\chi \approx 0.083$.</figcaption></figure>

### 3.2. 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 $\ell_{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*₂, 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*₁. 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 *l*₀ in Equation (16) for 2D models is given by

$$l_0 = 2.43a \ (2D \ rigid, \ XL5), \quad l_0 = 2.45a \ (2D \ empty, \ XL5). \qquad (17)$$

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*₁⟩, ⟨*S*₂⟩, and 3*N'* 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.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) The stress τsim vs. strain ε obtained by the 2D rigid model for the experimental data XL5 (×) and (b) the corresponding χ vs. ε; (c) τsim vs. ε obtained by the 2D empty model and (d) the corresponding χ vs. ε. In both (a) and (c), the symbol (○) denotes the data without crystallization obtained by assuming the same values of parameters *b* and *κ* assumed for the data (△).</figcaption></figure>
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ᵢ, as mentioned in Section 2.2.

In Figure 10a–f, $\chi$ 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 $\chi$ has an intermediate value between 0 and $\chi_{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 $\chi$ has an intermediate value between $\chi = 0$ and $\chi_{max}$ in Figure 10b,e. This observation indicates that $\chi$ smoothly changes with respect to strain $\varepsilon$ 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 $\sigma$ is performed only once every 100 MCSs. As mentioned in Section 2.4, the crystallization ratio $\chi$ at the large strain region becomes slightly smaller than $\chi_{max}$ if the MC update is performed more frequently. This result comes from the fact that the convergence speed of the variable $\vec{r}$ is very slow compared with that of $\sigma$. In other words, if $\sigma$ is updated more frequently before the lattice configuration is equilibrated, then the crystalline bonds of $\sigma = 1$ are expected to be distributed non-uniformly, and $\chi$ will be smaller than $\chi_{max}$ at the large strain region.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** The crystallization ratio $\chi$ vs. Monte Carlo sweep (MCS) of the 2D rigid model obtained at (a) $\varepsilon \approx 2.65$, (b) $\varepsilon \approx 3.47$, and (c) $\varepsilon \approx 7.58$, and the ratio $\chi$ vs. MCS of the 2D empty model at (d) $\varepsilon \approx 2.63$, (e) $\varepsilon \approx 3.45$, and (f) $\varepsilon \approx 5.90$.</figcaption></figure>

### 3.3. 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 viewpoint, 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 $\tau_{sim}$ and $\chi$ 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, \kappa$) are fixed
to $(b, \kappa) = (0.0015, 0.3)$ and $(b, \kappa) = (0.003, 0)$ for the experimental data XL5 and XL10, respectively. The value of $\ell_0$ in Equation (16) for the 3D model is given by

$$ \ell_0 = 2.69a \ (3D \ empty, \ XL5), \quad \ell_0 = 2.84a \ (3D \ empty, \ XL10). \qquad (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 $\chi$ 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 $\kappa$ is fixed to $\kappa = 0.85$. In contrast, the 3D lattice is composed of tetrahedrons, and, hence, the vertices hardly move compared to those on the 2D lattice.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) The stress $\tau_{sim}$ vs. strain $\varepsilon$ obtained by the 3D empty model for the experimental data XL5 (×) and (b) the corresponding $\chi$ vs. $\varepsilon$; (c) $\tau_{sim}$ vs. $\varepsilon$ obtained by the 3D empty model for the experimental data XL10 (×) and (d) the corresponding $\chi$ vs. $\varepsilon$. In both (a) and (c), the symbol (○) denotes the data without crystallization obtained by assuming the same values of parameters $b$ and $\kappa$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** The crystallization ratio $\chi$ vs. MCS of the 3D empty model obtained at (a) $\varepsilon \approx 1.37$, (b) $\varepsilon \approx 2.49$, and (c) $\varepsilon \approx 8.43$.</figcaption></figure>
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⁻¹⁰ [m]). From Table 1, we find that all assumed *a* values are larger than 10⁻¹⁰ [m], and, therefore, lattice modeling is meaningful as a coarse-grained technique for polymer networks.

<table><caption>**Table 1.** The lattice spacing *a* used for $\tau_{sim}$ (in Equation (14)) plotted in Figures 9 and 11.</caption><thead><tr><th>Figures</th><th>Model</th><th>a[m] for Data (○)</th><th>a[m] for Data (△)</th></tr></thead><tbody><tr><td>Figure 9a</td><td>2D rigid</td><td>1.305 × 10<sup>-9</sup></td><td>1.307 × 10<sup>-9</sup></td></tr><tr><td>Figure 9c</td><td>2D empty</td><td>1.329 × 10<sup>-9</sup></td><td>1.365 × 10<sup>-9</sup></td></tr><tr><td>Figure 11a</td><td>3D empty</td><td>1.265 × 10<sup>-9</sup></td><td>1.269 × 10<sup>-9</sup></td></tr><tr><td>Figure 11c</td><td>3D empty</td><td>1.654 × 10<sup>-9</sup></td><td>1.629 × 10<sup>-9</sup></td></tr></tbody></table>

Finally, in this subsection, we plot the mean length $\ell_{cr}$ for crystalline bonds in Figure 13a,b. This $\ell_{cr}$ is not always identical to $\ell_i^c$ in Equation (7) because of the constraint $\ell_{cr} < \ell_i^c$ for the crystalline bond length in the empty bond model. However, as described in Section 2.2, $\ell_i^c$ is dynamically changeable, and, hence, $\ell_{cr}$ is also dynamically changeable and depends on the strain $\epsilon$. The sudden jump of $\ell_i^c$ in Figure 13a for the 2D models implies that there are no crystalline bonds in the small strain region close to $\epsilon \approx 0$, and the nonzero $\ell_i^c$ 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 $\ell_i^c$ increases with increasing $\epsilon$, and this behavior is consistent with that observed in real materials [3].

<figure><img src="image_15.png" /><figcaption>**Figure 13.** The mean length $\ell_{cr}$ of crystalline bonds corresponding to (a) 2D rigid and empty models and (b) the 3D empty model for the experimental data XL5 and XL10.</figcaption></figure>

### 3.4. Dependence on Stiffness $\kappa$

In the preceding subsection, two different experimental data, XL5 and XL10, are simulated and plotted in Figure 11a,c, where the stiffness $\kappa = 0.3$ and $\kappa = 0$ are assumed, respectively. However, the dependence of the results on $\kappa$ is not so clear. Therefore, in this subsection, we should discuss the effect of stiffness $\kappa$ on the stress–strain curve and clarify that the effect of $\kappa$ is different from the effect of SIC. First, we plot in Figure 14a the results of 3D model with the crystalline state (○) ($\Leftrightarrow \chi > 0$) and those without the crystalline state (△, ▽) ($\Leftrightarrow \chi = 0$), where $\kappa$ is fixed to $\kappa = 1$ and $\kappa = 0$, respectively. In Figure 14a, the values of $\tau$ 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 $\chi = 0$, except in the smaller strain region. This decrement of $\tau$ (○) is understood to be due to the crystalline states.

To see the effect of $\kappa$ more clearly, the stress $\tau_{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 $\kappa$ increases the stress only at the region $\epsilon \le 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 $\varepsilon \ge 2.5$, for example. Thus, we find from these numerical data that the stress relaxation by SIC is independent of the effect of $\kappa$. 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 $\kappa$ 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) and (5).

<figure><img src="image_16.png" /><figcaption>**Figure 14.** (a) The stress $\tau$ vs. $\varepsilon$ obtained by the 3D model with the crystalline state (○) and without the crystalline state (Δ, ∇), and (b) the stress $\tau_{\text{sim}}$ vs. $\varepsilon$ obtained without the crystalline state. The stress $\tau$ 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 $\tau_{\text{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.</figcaption></figure>

**4. 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 $\sigma$ 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 $\chi$ 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.
**Author Contributions:** H.K. designed the study, performed the modeling and simulations, and wrote the paper. V.E. performed the simulations. C.B., J.-M.C., G.D., G.S., L.C., J.-Y.C., and T.T. evaluated the simulation results and data analyses and critically reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported in part by the Collaborative Research Project of the Institute of Fluid Science, Tohoku University, and JSPS KAKENHI Grant No. JP17K05149.

**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

## References

1. Toki, S.; Sics, I.; Hsiao, B.S.; Murakami, S.; Tosaka, M.; Poompradub, S.; Kohjiya, S.; Ikeda, Y. Structural developments in synthetic rubbers during uniaxial deformation by in situ synchrotron X-ray diffraction. *J. Polym. Sci. B* **2004**, *42*, 956–964. [CrossRef]

2. Albouy, P.-A. The conformation of poly (dimethylsiloxane) in the crystalline state. *Polymer* **2000**, *41*, 3083–3086. [CrossRef]

3. Albouy, P.-A.; Sotta P. Strain-Induced Crystallization in Natural Rubber. In *Polymer Crystallization II*; Auriemma, F., Alfonso, G.C., Rosa C.D., Eds.; Springer: Cham, Switzerland, 2017, pp. 167-205.

4. Imbernon, L.; Puachet, R.; Pire, M.; Albouy, P.-A.; Tencé-Giraut, S.; Norves, S. Strain-induced crystallization in sustainably crosslinked epoxidized natural rubber. *Polymer* **2016**, *93*, 189–197. [CrossRef]

5. Candau N.; Laghmach, R.; Chazeau, L.; Chenal, J.M.; Gauthier, C.; Biben, T.; Munch, E. Influence of strain rate and temperature on the onset of strain induced crystallization in natural rubber *Eur. Polym. J.* **2015**, *64*, 244–252.

6. Urayama, K. An experimentalist's view of the physics of rubber elasticity. *J. Polym. Sci. B* **2006**, *44*, 3440–3444. [CrossRef]

7. Kuhn, W. Beziehungen zwischen Molekülgröse, statistischer Molekülgestalt und elastischen Eigenschaften hochpolymerer Stoffe. *Kolloid Z.* **1936**, *76*, 258–271. [CrossRef]

8. Kuhn, W. Dependence of the average transversal on the longitudinal dimensions of statistical coils formed by chain molecules. *J. Polym. Sci.* **1946**, *1*, 380–388. [CrossRef]

9. Treloar, L.R.G. Rubber elasticity. *Contemp. Phys.* **1971**, *12*, 33–56. [CrossRef]

10. Treloar, L.R.G. *The Physics of Rubber Elasticity*, 3rd ed.; Clarendon Press: Oxford, UK, 1975.

11. Flory, P.J. Network topology and the theory of rubber elasticity. *Polym. Int.* **1985**, *17*, 96–102. [CrossRef]

12. Tanaka, F. *Polymer Physics: Applications to Molecular Association and Thermoreversible Gelation*; Cambridge University Press Ithaca: Cambridge, UK, 2011.

13. Doi, M.; Edwards, S.F. *The Theory of Polymer Dynamics*; Oxford University Press: Oxford, UK, 1986.

14. Bowick, M.; Travesset, A. The statistical mechanics of membranes. *Phys. Rep.* **2001**, *344*, 255–308. [CrossRef]

15. Wiese, K.J. Polymerized Membranes, a Review. In *Phase Transitions and Critical Phenomena* **19**; Domb, C., Lebowitz, J.L., Eds.; Academic Press: London, UK, 2000; pp. 253–498.

16. Nelson, D. The Statistical Mechanics of Membranes and Interfaces. In *Statistical Mechanics of Membranes and Surfaces*, 2nd ed.; Nelson, D., Piran, T., Weinberg, S., Eds.; World Scientific: Singapore, 2004; pp. 1–17.

17. Gompper, G.; Kroll, D.M. Triangulated-surface Models of Fluctuating Membranes. In *Statistical Mechanics of Membranes and Surfaces*, 2nd ed.; Nelson, D., Piran, T., Weinberg S., Eds.; World Scientific: Singapore, 2004; pp. 359–426.

18. Deam, R.T.; Edwards, F.R.S. The theory of rubber elasticity. *Phil. Trans. R. Soc. Lond.* **1976**, *280*, 317–353.

19. Immirzi A.; Tedesco C.; Monaco G.; Tonelli, A.E. Crystal structure and melting entropy of natural rubber. *Macromolecules* **2005**, *38*, 1223–1231. [CrossRef]
20. Koibuchi, H.; Bernard, C.; Chenal, J.M.; Diguet, G.; Sebald, G.; Cavaille, J.Y.; Takagi, T.; Chazeau, C. Monte Carlo Study of Rubber Elasticity on the Basis of Finsler Geometry Modeling. *Symmetry* **2020**, *11*, 1124. [CrossRef]

21. Metropolis, N.; Rosenbluth, A.W.; Rosenbluth, M.N.; Teller, A.H. Equation of State Calculations by Fast Computing Machines. *J. Chem. Phys.* **1953**, *21*, 1087–1092. [CrossRef]

22. Landau, D.P. Finite-size behavior of the Ising square lattice. *Phys. Rev. B* **1976**, *13*, 2997–3011. [CrossRef]

23. Wheater J.F. Random surfaces: From polymer membranes to strings. *J. Phys. A Math. Gen.* **1994**, *27*, 3323–3354. [CrossRef]

24. Pradhan, S.; Costa, F.R.; Wagenknecht, U.; Jehnichen, D.; Bhowmick, A.K.; Heinrich, G. Elastomer/LDH nanocomposites: Synthesis and studies on nanoparticle dispersion, mechanical properties and interfacial adhesion. *Eur. Polym. J.* **2008**, *44*, 3122–3132. [CrossRef]

25. Mooney, M. A Theory of Large Elastic Deformation. *J. Appl. Phys.* **1940**, *11*, 582–592. [CrossRef]

26. Rivlin, R.S. Large elastic deformations of isotropic materials I. Fundamental concepts. *Phil. Trans. R. Soc. Lond. Soc. A* **1948**, *240*, 459–490.

<figure><img src="image_17.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).