Comparison of Linear Viscoelastic Behaviors of Langevin- and Jump-Rouse Models

Dynamics of unentangled polymer melts in or near equilibrium can be well described by the Rouse model[1, 2]. In the Rouse model, a single tagged polymer chain in the system is considered, and the tagged chain is expressed as an array of beads linearly connected by harmonic springs. Then the dynamics of the tagged chain is modeled by a simple Langevin equation (the Langevin-Rouse model). Although the Rouse model is originally developed as a model for dilute polymer solutions, it rather describes the dynamics of unentangled polymer melts. This is because the interaction between segments is screened[2, 3] and the hydrodynamic interaction is negligible in polymer melts.

The application of the Rouse model is not limited to simple unentangled polymer melts. For example, it can be used to describe the relaxation mode of an entangled polymers by the constraint-release (CR) mechanism[4, 5]. In the simple reptation model[2], a single tagged entangled polymer chain is modeled as a chain confined in a tube-like constraint, and the tagged chain can relax by the reptation motion along the tube and by the contour length fluctuation. Because the tube is formed by the surrounding chains, the relaxation of the surrounding chains induce some local motions of the tube. Such local motions may be interpreted as the local rearrangement of tube segments described by the Rouse-like dynamics (the CR-Rouse model)[4]. In the CR-Rouse model, the tube segment moves by a local jump[6], not by the Langevin equation.

The dynamics of unentangled polymers with associative groups is sometimes described by a similar Rouse-like model[7, 8, 9]. Associative groups effectively work as transient cross-links. An associative group in a cross-link can dissociate from the cross-link by the thermal fluctuation, and then reassociate to another cross-link. The motion of an associative group may look as discrete jumps. Then the chain dynamics can be modeled by the Rouse-like dynamics (the sticky Rouse model). The dynamics of unentangled supercooled polymers[10, 11] may also be (roughly) described by a similar model. In a supercooled liquid, segmental motion becomes highly correlated and constrained (the dynamic heterogeneity)[12]. If we observe the motion of a single particle in a supercooled liquid, at the short-time scale it is effective trapped in a cage formed by surrounding particles. At the long-time scale, the particle escapes from the cage and then trapped to another cage again[12, 13, 14]. At the long-time scale, the segmental motion may be modeled by discrete jumps rather than continuum Brownian motion. If we apply this picture to a supercooled unentangled polymer chain, we will have the Rouse-like dynamics driven by local jumps.

Naively, we expect that the chain dynamics of such a jump-based model (the jump-Rouse model) is similar to that of the Langevin-Rouse model. Intuitively, the linear viscoelasticity of the jump-Rouse model will be approximated well by that of the Langevin-Rouse model. However, because these two models are not equivalent, there is no guarantee that the dynamic properties of two models are the same even for the same parameter set. Direct comparison of linear viscoelasticity data by two models with the same parameter set will be required to judge whether the linear viscoelasticity of the jump-Rouse model can be approximated by the Langevin-Roue model or not.

The purpose of this work is to investigate the linear viscoelasticity of the Rouse model described by the jump-Rouse model in detail. To compare the Langevin- and jump-Rouse models, we develop a jump rate model with a freely tunable resampling ratio. Our model reduces to the Langevin-Rouse model at a certain limit, and this property allows us to connect two models smoothly and compare dynamic properties of two models directly. We perform numerical simulations for the jump-Rouse model and calculate the shear relaxation modulus. Then we compare the shear relaxation moduli of the jump-Rouse model with those of the Langevin-Rouse model. We show that the jump-Rouse model exhibits similar but slightly different linear viscoelasticity, compared with that of the Langevin-Rouse model. The jump-Rouse model exhibits a characteristic relaxation at the short-time scale. At the long-time scale, the Langevin- and jump-Rouse models give almost the same shear relaxation modulus. We discuss the origin of the short-time relaxation mode in the jump-Rouse model. Our result justifies the use of the shear relaxation moduli by the Langevin-Rouse model as an approximate expression for that by the jump-Rouse model.

We perform simulations for the jump-Rouse model, in order to calculate the linear viscoelasticity. In this work, we simulate the jump process by using the kinetic Monte Carlo (kMC) method[17, 18]. The kMC method enables us to simulate the stochastic process exactly.

In what follows, we use the dimensionless units. We set $k_{B}T=1$, $b=1$, and $3k_{B}T/\zeta b^{2}=1$. Also, we set the bead density to unity, $\rho=1$, because the bead density is not that important in our model. In what follows, we express all the physical quantities in the dimensionless units.

In the kMC method, the time-evolution of the bead positions $\{\bm{R}_{i}(t)\}$ is described by sequence of resampling events. We describe the time where the $i$-th resampling event occurs as $t_{i}$ ($i=1,2,3,\dots$). For convenience, we set $t_{0}=0$. The resampling events occur stochastically, and thus the time step between two successive events $t_{i}-t_{i-1}$ is a stochastic quantity. $t_{i}-t_{i-1}$ obeys the exponential distribution with $\langle t_{i}-t_{i-1}\rangle=1/Nk_{\text{jump}}$. In dimensionless units, $k_{\text{jump}}=2/[1-(1-\theta)^{2}]$, and thus $t_{i}$ is sampled as:

where $\epsilon$ is a random number sampled from the exponential distribution with the first order moment $\langle\epsilon\rangle=1$. At $t=t_{i}$, the target bead index is randomly selected from $j=1,2,\dots,N$, with the equal probability. Then the position of the selected bead is resampled. If we assume that the $j$-th bead is selected, its position is resampled as

with $\alpha_{j}$ and $\beta_{j}$ defined by eqs (21) and (22). $\bm{w}$ is a random number sampled from the Gaussian random number of which moments are $\langle\bm{w}\rangle=0$ and $\langle\bm{w}\bm{w}\rangle=\bm{1}$. Other bead positions are unchanged by this resampling: $\bm{R}_{k}(t_{i+1})=\bm{R}_{k}(t_{i})$ ($k\neq j$). By iterating the resampling events, we can generate the time-dependent bead position $\{\bm{R}_{j}(t)\}$.

The kMC scheme explained above is implemented on Octave (version 5.2.0)[19], and the simulations are performed with various parameters. We vary the number of beads per chain $N$ as $N=5,10,20,40,$ and $80$. We also vary the resampling ratio $\theta$ as $\theta=0.0625,0.125,0.25,0.5,$ and $1$. For all the cases, simulations with $128$ different random seeds are performed. Each simulation is run upto $t=10^{5}$. The initial bead positions $\{\bm{R}_{j}(0)\}$ are sampled from the equilibrium distribution.

During the simulation, the time series data of the stress tensor are stored, and the relaxation moduli are calculated from the stress tensor data. The single-chain stress tensor at time $t$ is calculated as

The linear response theory[20, 21] gives the expression of the shear relaxation modulus in terms of the shear stress correlation:

where and $\langle\dots\rangle_{\text{eq}}$ is the equilibrium ensemble average. To improve the statistical accuracy, we use the following formula by Likhtman[21] instead of eq (30):

where $\hat{N}_{\alpha\beta}(t)=\hat{\sigma}_{\alpha\alpha}(t)-\hat{\sigma}_{\beta\beta}(t)$ is the normal stress difference. We calculate the ensemble average in eq (31) as the averages over time and random seeds.

Figure 2 shows the shear relaxation moduli of the jump-Rouse model with various $N$ and $\theta$. We observe that the shear relaxation modulus is not sensitive to $\theta$. We estimate the viscoelastic Rouse time $\tau_{R}$ from the shear relaxation modulus at the long-time region:

The viscoelastic Rouse times estimated by the least-squares fitting are shown in Figure 3. For comparison, the viscoelastic Rouse time for the Langevin-Rouse model, $\tau_{R}=1/8\sin^{2}(\pi/2N)$ is also shown in Figure 3. As expected from Figure 2, the Rouse time of the jump-Rouse model is not sensitive to $\theta$, and almost coincides to that of the Langevin-Rouse model. We expect that the shear relaxation modulus by the jump-Rouse model can be well approximated by that of the Langevin-Rouse model.

Before we proceed to the detailed analyses for the relaxation modulus data by the jump-Rouse model, here we discuss the terminal region. If $N$ is sufficiently large, the Langevin-Rouse model exhibits universal relaxation behavior which is independent of $N$. This property makes the Langevin-Rouse model physically reasonable. To examine whether the jump-Rouse model exhibits the universal relaxation behavior at the terminal region or not, we reduce $t$ and $G(t)$ as $t/\tau_{R}$ and $NG(t)$ and compare the reduced data. Figure 4 shows the reduced relaxation modulus data for $N=40$ and $80$, and $\theta=0.0625,0.25,$ and $1$. We observe that all the data are almost collapsed. Moreover, we observe that the power-law type relaxation behavior, $G(t)\propto N^{-1/2}$, the behavior well-known for the Langevin-Rouse model for sufficiently large $N$. Thus we find that the jump-Rouse model exhibits the same universal behavior as the jump-Rouse model, and it can be employed as a physically reasonable model in the same way as the Langevin-Rouse model.

If we observe the $\theta$-dependence in Figure 2 in detail, we find that the shear relaxation modulus changes systematically by changing $\theta$. As increasing $\theta$, a weak shoulder appears at the short-time region ($t\lesssim 10^{0}$). In addition, for small $N$ cases, the terminal relaxation seems to be slightly accelerated. At the limit of $\theta\to 0$, the jump-Rouse model should reduce to the Langevin-Rouse model, as we explained. Thus we expect that the data for $\theta=0.0625$ are close to the shear relaxation modulus by the Langevin-Rouse model. Also, we expect that the jump-Rouse model with $\theta=1$ is the most deviated from the Langevin-Rouse model (although the deviation may be rather small). Thus we compare the shear relaxation moduli for $\theta=0.0625$ and $1$, with those by the Langevin-Rouse model. Figure 5 shows the shear relaxation moduli for $\theta=0.0625$ and $1$, taken from Figure 2, and the shear relaxation modulus of the Langevin-Rouse model. In dimensionless units, the shear relaxation modulus by eq (4) becomes as follows:

We find that the data for $\theta=0.0625$ are almost perfectly overlap to eq (33). Thus we confirm that our jump-Rouse model reduces to the Langevin-Rouse model if $\theta$ is sufficiently small.

The data for $\theta=1$ is deviated from the Langevin-Rouse model, and the deviation is large at the short-time region ($t\lesssim 10^{0}$). The jump-Rouse model exhibits a shoulder whereas the Langevin-Rouse model exhibits a power-law type relaxation ($G(t)\propto t^{-1/2}$). The shoulder at the short-time region can be interpreted as the relaxation mode due to the local jumps. In the case of $\theta=1$, the local bead position is fully relaxed by a single jump event. Also, if the chain is sufficiently long ($N\gg 1$), the contribution of the chain ends are negligibly small. Therefore, we expect that we have a short-time relaxation mode which is almost independent of $N$. The local jumps are not correlated, and thus there is no collective mode at the short-time scale. This is why the jump-Rouse model does not exhibit the power-law type relaxation, which originates from the higher order Rouse modes.

We estimate the short-time shear relaxation modulus based on jump events. We describe the initial bead positions as $\{\bm{R}_{j}\}$. We consider the case where the $j$-th bead is resampled by the first jump. We assume that $N$ is sufficiently large, and assume that $2\leq j\leq N-1$. The single-chain stress at $t=0$ can be decomposed into the $\bm{R}_{j}$-independent and $\bm{R}_{j}$-dependent parts:

Here, $\hat{\bm{\sigma}}^{\prime}_{j}(\{\bm{R}_{j}\})$ is the $\bm{R}_{j}$-independent part of the stress tensor:

The stress tensor is unchanged until the first jump event, and thus $\hat{\bm{\sigma}}(t)=\hat{\bm{\sigma}}(0)$ for $t<t_{1}$. At the first jump, only the $\bm{R}_{j}$-dependent part of the stress is changed. If we describe the newly sampled bead position as $\bm{R}_{j}^{\prime}$, the stress tensor at $t_{1}\leq t<t_{2}$ becomes

If $t$ is sufficiently small, we will observe at most one jump event from time $0$ to time $t$. If there is no jump event ($t<t_{1}$), the correlation of the shear stress is simply calculated as

If there is just one jump event ($t_{1}\leq t<t_{2}$), the correlation of the shear stress can be still analytically calculated. In equilibrium, $j$-dependent and $j$-independent parts of the stress tensor are not correlated each other. In addition, the $j$-dependent part can be evaluated by using the local equilibrium distribution (9):

The probability that the $j$-th bead position is not resampled upto time $t$ is estimated as $\exp(-k_{\text{jump}}t)=\exp(-2t)$. Also, the probability to have two or more jump events will be sufficiently small and negligible. Thus the average stress correlation becomes

To derive eq (39), we considered a jump event only for the $j$-th bead. There are $N$ beads in the chain, and all the beads can jump with an equal probability. Jump of another bead reduces the stress correlation in the same way as eq (39). By considering the jump events for all the beads, the stress correlation is estimated as

Finally we have the following simple expression for the shear relaxation modulus:

Here we have assumed that $N$ is sufficiently large and used the approximation $(N-1)/N\approx 1$. Eq (41) means that the shear relaxation modulus of the jump-Rouse model with $\theta=1$ is well approximated by the single Maxwell form. We observe that eq (41) reasonably explain the simulation data in Figure 5(b). As we discussed, the jump-Rouse model exhibits universal terminal relaxation behavior if $N$ is sufficiently large. After the short-time relaxation mode discussed here is relaxed, we will simply observe the relaxation modulus which is almost the same as that by the Langevin-Rouse model.

The characteristic short-time relaxation behavior by eq (41) will be, however, not important in most of practical cases. In the case of entangled polymers, there are various relaxation modes such as the reptation and the contour length fluctuation modes. The associative polymers have short-time subchain relaxation. In the case of supercooled unentangled polymers, we have the glassy relaxation mode at the short-time region. Then, we will observe other relaxation modes at the time scale where the shoulder is expected to emerge. A weak shoulder will be apparently mixed with other relaxation modes, and then we will not observe a clear shoulder. At the long-time region, the jump-Rouse model gives approximately the same shear relaxation modulus as the Langevin-Rouse model. Therefore, we conclude that our simulation result justifies the use of the shear relaxation modulus of the Langevin-Rouse as a good approximation for that of the jump-Rouse model.

Before we end this section, we briefly discuss the temperature dependence of the jump-Rouse model. If the temperature of the system is changed, the relaxation time and the characteristic modulus in the Langevin-Rouse model are changed but the relaxation mode distribution itself is not changed. Thus the time-temperature superposition (tTS) holds for the Langevin-Rouse model. In contrast to the Langevin-Rouse model, the temperature dependence of the jump-Rouse model is not that clear. The thermal activation process in the jump-Rouse model is expressed by the resampling ratio $\theta$, and how it depends on the temperature is not clear (while the temperature dependence of the thermal noise $\bm{\xi}_{j}(t)$ for the Langevin-Rouse model is rather simple). It would be natural to expect that $\theta$ decreases as the temperature increases. Then the shape of the relaxation modulus at the short-time region is expected to change as the temperature changes. This means that the tTS will not hold for the jump-Rouse model. We may be able to extract some information on the local jumps from the failure of the tTS. But as we discussed, the shoulder in the short-range region will not be clearly observed in practice. If we observe only on the long-time region, the jump-Rouse model exhibits the universal behavior which is independent of $\theta$. Thus we expect that the tTS will approximately hold for the jump-Rouse model, in most of practical cases.

We constructed a jump rate model with a variable resampling ratio for the jump-Rouse model, and performed numerical simulations. By tuning the resampling ratio, the jump-Rouse model reduces both to the Langevin-Rouse model and the local resampling model. We calculated the shear relaxation modulus data from by simulations, and investigated the effects of the number of beads $N$ and the resampling ratio $\theta$. The effect of $\theta$ on the relaxation modulus is weak. For small $\theta$ (such as $\theta=0.0625$), the shear relaxation modulus by the jump-Rouse model almost coincides to that by the Langevin-Rouse model. For $\theta=1$, we observe small deviations. At the short-time region where the Langevin-Rouse model exhibits the power-law type relaxation, we observe a small shoulder in the jump-Rouse model. Also, if $N$ is small, the viscoelastic relaxation time is slightly accelerated in the jump-Rouse model. We theoretically analyzed the short-time relaxation mode in the jump-Rouse model with $\theta=1$, and derived that the relaxation modulus becomes a single Maxwell model for $N\gg 1$. Except the short-time region, the shear relaxation modulus for $\theta=1$ is close to that by the Langevin-Rouse model. We consider that our result justifies the use of the shear relaxation modulus of the Langevin-Rouse model as a good approximation for that of the jump-Rouse model.

This work was supported by JST, PRESTO Grant Number JPMJPR1992, Japan, Grant-in-Aid (KAKENHI) for Scientific Research Grant B No. JP19H01861, and Grant-in-Aid (KAKENHI) for Scientific Research Grant B No. JP23H01142.