A Transient Bond Model for Dynamic Constraints in Meso-Scale Coarse-Grained Systems

TL;DR

This paper introduces a highly coarse-grained transient bond model for simulating the dynamic constraints in entangled polymers, effectively reproducing their viscoelastic and diffusion behaviors with simplified particles and bonds.

Contribution

The paper develops a novel coarse-grained slip-spring model based on transient bonds that captures the dynamics of entangled polymers without affecting static properties.

Findings

Model accurately reproduces viscoelastic behavior of entangled polymers

Transient bonds influence dynamics but not static equilibrium properties

Simulation results match experimental observations of polymer diffusion

Abstract

The dynamical properties of entangled polymers originate from the dynamic constraints due to the uncrossability between polymer chains. We propose a highly coarse-grained simulation model with transient bonds for such dynamically constrained systems. Based on the ideas of the responsive particle dynamics (RaPiD) model [P. Kindt and W. J. Briels, J. Chem. Phys. 127, 134901 (2007)] and the multi-chain slip-spring model [T. Uneyama and Y. Masubuchi, J. Chem. Phys. 137, 154902 (2012)], we construct the RaPiD type transient bond model as a coarse-grained slip-spring model. In our model, a polymer chain is expressed as a single particle, and particles are connected by transient bonds. The transient bonds modulate the dynamics of particles but they do not affect static properties in equilibrium. We show the relation between parameters for the entangled polymer systems and those for the transient bond model. By performing simulations based on the transient bond model, we show how model parameters affect the linear viscoelastic behavior and the diffusion behavior. We also show that the viscoelastic behavior of entangled polymer systems can be well reproduced by the transient bond model.

Equations

$$\mathcal{Z}=\frac{1}{{\Lambda }^{3M}M!}\int d\{{\mathbf{R}}_i\}\exp \left[-\frac{1}{k_{B}T}\sum _{i>j}v({\mathbf{R}}_i-{\mathbf{R}}_j)\right],$$

$$P_{\text{eq}}(\{{\mathbf{R}}_i\})=\frac{1}{\mathcal{Z}{\Lambda }^{3M}M!}\exp \left[-\frac{1}{k_{B}T}\sum _{i>j}v({\mathbf{R}}_i-{\mathbf{R}}_j)\right].$$

$${\mathcal{U}}_{\text{transient}}(\{{\mathbf{R}}_i\},\{n_{ij}\})=\sum _{i>j}{n}_{ij}u({\mathbf{R}}_i-{\mathbf{R}}_j).$$

$$\sum _{\{n_{ij}\}}{P}_{\text{eq}}(\{{\mathbf{R}}_i\},\{n_{ij}\})=P_{\text{eq}}(\{{\mathbf{R}}_i\}),$$

$$P_{\text{eq}}(\{n_{ij}\}|\{{\mathbf{R}}_i\})=\frac{1}{\Xi (\{{\mathbf{R}}_i\})\prod _{i>j}{n}_{ij}!}\exp \left[\frac{\mu }{k_{B}T}\sum _{i>j}{n}_{ij}-\frac{1}{k_{B}T}\sum _{i>j}{n}_{ij}u({\mathbf{R}}_i-{\mathbf{R}}_j)\right].$$

$$\begin{array}{r}\begin{array}{rl}\Xi (\{{\mathbf{R}}_i\})& =\prod _{i>j}\sum _{n_{ij}=0}^{\infty }\frac{1}{n_{ij}!}\exp \left[\frac{\mu }{k_{B}T}{n}_{ij}-\frac{1}{k_{B}T}{n}_{ij}u({\mathbf{R}}_i-{\mathbf{R}}_j)\right]\\ & =\exp \left[\xi \sum _{i>j}{e}^{-u({\mathbf{R}}_i-{\mathbf{R}}_j)/k_{B}T}\right],\end{array}\end{array}$$

$$\begin{split}P_{\text{eq}}(\{\bm{R}_{i}\},\{n_{ij}\})&=P_{\text{eq}}(\{n_{ij}\}|\{\bm{R}_{i}\})P_{\text{eq}}(\{\bm{R}_{i}\})\\ &=\frac{1}{\displaystyle\mathcal{Z}\Lambda^{3M}M!}\left[{\prod_{i>j}\frac{\xi^{n_{ij}}}{n_{ij}!}}\right]\exp\Bigg{[}-\frac{1}{k_{B}T}\sum_{i>j}n_{ij}u(\bm{R}_{i}-\bm{R}_{j})\\ &\qquad-\frac{1}{k_{B}T}\sum_{i>j}v(\bm{R}_{i}-\bm{R}_{j})-{\xi}\sum_{i>j}e^{-u(\bm{R}_{i}-\bm{R}_{j})/k_{B}T}\Bigg{]}.\end{split}$$

$$\mathcal{F}(\{{\mathbf{R}}_i\},\{n_{ij}\})\equiv \sum _{i>j}v({\mathbf{R}}_i-{\mathbf{R}}_j)+\sum _{i>j}{n}_{ij}u({\mathbf{R}}_i-{\mathbf{R}}_j)+k_{B}T\xi \sum _{i>j}{e}^{-u({\mathbf{R}}_i-{\mathbf{R}}_j)/k_{B}T}.$$

$$\frac{d{\mathbf{R}}_i(t)}{dt}=-\frac{1}{\zeta }\frac{\partial \mathcal{F}(\{{\mathbf{R}}_i\},\{n_{ij}\})}{\partial {\mathbf{R}}_i}+\mathbf{\kappa }(t)\cdot {\mathbf{R}}_i+\sqrt{\frac{2k_{B}T}{\zeta }}{\mathbf{w}}_i(t),$$

$$⟨{\mathbf{w}}_i(t)⟩=0,⟨{\mathbf{w}}_i(t){\mathbf{w}}_j(t^{\prime })⟩={\delta }_{ij}\delta (t-t^{\prime })1,$$

$$W(n_{ij}-1|n_{ij})=\frac{n_{ij}}{\tau },$$

$$W(n_{ij}|n_{ij}+1)P_{\text{eq}}(n_{ij}+1|\{{\mathbf{R}}_i\})=W(n_{ij}+1|n_{ij})P_{\text{eq}}(n_{ij}|\{{\mathbf{R}}_i\}).$$

$$W(n_{ij}+1|n_{ij})=\frac{1}{\tau }\xi e^{-u({\mathbf{R}}_i-{\mathbf{R}}_j)/k_{B}T}.$$

$$\begin{array}{rl}\hat{\mathbf{\sigma }}& =\frac{1}{V}\sum _{i>j}\frac{\partial \mathcal{F}(\{{\mathbf{R}}_i\},\{n_{ij}\})}{\partial ({\mathbf{R}}_i-{\mathbf{R}}_j)}({\mathbf{R}}_i-{\mathbf{R}}_j)-\frac{M{k}_{B}T}{V}1\\ & =\frac{1}{V}\sum _{i>j}\left[\left[n_{ij}-\xi e^{-u({\mathbf{R}}_i-{\mathbf{R}}_j)/k_{B}T}\right]\frac{\partial u({\mathbf{R}}_i-{\mathbf{R}}_j)}{\partial ({\mathbf{R}}_i-{\mathbf{R}}_j)}+\frac{\partial v({\mathbf{R}}_i-{\mathbf{R}}_j)}{\partial ({\mathbf{R}}_i-{\mathbf{R}}_j)}\right]({\mathbf{R}}_i-{\mathbf{R}}_j)-\frac{M{k}_{B}T}{V}1.\end{array}$$

$$\frac{\partial P(\{{\mathbf{R}}_i\},\{n_{ij}\},t)}{\partial t}=[{\mathcal{L}}_0+{\mathcal{L}}_1(t)]P(\{{\mathbf{R}}_i\},\{n_{ij}\},t).$$

$${\mathcal{L}}_1(t)P=-\sum _i\frac{\partial }{\partial {\mathbf{R}}_i}\cdot [\mathbf{\kappa }(t)\cdot {\mathbf{R}}_{i}P].$$

$$P(\{{\mathbf{R}}_i\},\{n_{ij}\},t)=P_{\text{eq}}(\{{\mathbf{R}}_i\},\{n_{ij}\})+P_1(\{{\mathbf{R}}_i\},\{n_{ij}\},t).$$

$$\begin{array}{rl}\frac{\partial P_1(\{{\mathbf{R}}_i\},\{n_{ij}\},t)}{\partial t}& \approx {\mathcal{L}}_1(t)P_{\text{eq}}(\{{\mathbf{R}}_i\},\{n_{ij}\})+{\mathcal{L}}_0{P}_1(\{{\mathbf{R}}_i\},\{n_{ij}\},t)\\ & =\frac{V}{k_{B}T}\hat{\mathbf{\sigma }}:\mathbf{\kappa }(t)P_{\text{eq}}(\{{\mathbf{R}}_i\},\{n_{ij}\})+{\mathcal{L}}_0{P}_1(\{{\mathbf{R}}_i\},\{n_{ij}\},t).\end{array}$$

$$P_1(\{{\mathbf{R}}_i\},\{n_{ij}\},t)=\frac{V}{k_{B}T}{\int }_{-\infty }^{t}d{t}^{\prime }{e}^{(t-t^{\prime }){\mathcal{L}}_0}\hat{\mathbf{\sigma }}:\mathbf{\kappa }(t^{\prime })P_{\text{eq}}(\{{\mathbf{R}}_i\},\{n_{ij}\}),$$

$$\begin{array}{rl}\mathbf{\sigma }(t)& =\int d\{{\mathbf{R}}_i\}\sum _{\{n_{ij}\}}\hat{\mathbf{\sigma }}[P_{\text{eq}}(\{{\mathbf{R}}_i\},\{n_{ij}\})+P_1(\{{\mathbf{R}}_i\},\{n_{ij}\},t)]\\ & =⟨\hat{\mathbf{\sigma }}{⟩}_{\text{eq}}+\frac{V}{k_{B}T}{\int }_{-\infty }^{t}d{t}^{\prime }{⟨\hat{\mathbf{\sigma }}(t-t^{\prime })\hat{\mathbf{\sigma }}⟩}_{\text{eq}}:\mathbf{\kappa }(t^{\prime }),\end{array}$$

$$G(t)=\frac{V}{k_{B}T}⟨{\hat{\sigma }}_{xy}(t){\hat{\sigma }}_{xy}{⟩}_{\text{eq}}.$$

$$v(\mathbf{r})$$

$$u(\mathbf{r})$$

$$\mathcal{Z}=\frac{V^M}{{\Lambda }^{3M}M!},$$

$$\frac{\mathcal{F}(\{{\mathbf{R}}_i\},\{n_{ij}\})}{k_{B}T}=\frac{3}{2a^2}\sum _{i>j}{n}_{ij}({\mathbf{R}}_i-{\mathbf{R}}_j{)}^2+\xi \sum _{i>j}{e}^{-3({\mathbf{R}}_i-{\mathbf{R}}_j{)}^2/2a^2}.$$

$$\hat{\mathbf{\sigma }}=\frac{3k_{B}T}{2a^{2}V}\sum _{i>j}\left[n_{ij}-\xi e^{-3({\mathbf{R}}_i-{\mathbf{R}}_j{)}^2/2a^2}\right]({\mathbf{R}}_i-{\mathbf{R}}_j)({\mathbf{R}}_i-{\mathbf{R}}_j)-\frac{M{k}_{B}T}{V}1.$$

$$⟨\hat{\mathbf{\sigma }}{⟩}_{\text{eq}}=\int d\{{\mathbf{R}}_i\}\sum _{\{n_{ij}\}}\hat{\mathbf{\sigma }}{P}_{\text{eq}}(\{n_{ij}\}|\{{\mathbf{R}}_i\})P_{\text{eq}}(\{{\mathbf{R}}_i\})=-\frac{M{k}_{B}T}{V}1.$$

$$\begin{array}{rl}⟨n_{ij}{⟩}_{\text{eq}}& =\int d{\mathbf{R}}_{i}d{\mathbf{R}}_j\sum _{n_{ij}}{n}_{ij}{P}_{\text{eq}}({\mathbf{R}}_i,{\mathbf{R}}_j,n_{ij})\\ & =\int d{\mathbf{R}}_{i}d{\mathbf{R}}_j\frac{1}{V^2}\sum _{n_{ij}=0}^{\infty }\frac{n_{ij}}{n_{ij}!}{\xi }^{n_{ij}}\exp \left[-\frac{3n_{ij}({\mathbf{R}}_i-{\mathbf{R}}_j{)}^2}{2a^2}-\xi e^{-3({\mathbf{R}}_i-{\mathbf{R}}_j{)}^2/2a^2}\right]\\ & =\frac{\xi }{V}{\left(\frac{2\pi a^2}{3}\right)}^{3/2}.\end{array}$$

$$\varphi =\frac{M(M-1)}{2V}⟨n_{ij}{⟩}_{\text{eq}}\approx \frac{\xi {\rho }^2}{2}{\left(\frac{2\pi a^2}{3}\right)}^{3/2},$$

$$\rho =\frac{{\rho }_{\text{segment}}}{N}.$$

$$a^2\propto N{b}^2.$$

Full text

I Introduction

The entangled polymer systems exhibit characteristic relaxation behavior such as very long relaxation time which strongly depends on the degree of polymerizationDoi and Edwards (1986). Because of their very long relaxation times, the simulations for polymer melts and solutions with large degrees of polymerization are difficult. Especially, simulations for long time relaxation relaxation processes by microscopic molecular dynamics models (such as the Kremer-Grest molecular dynamics modelKremer and Grest (1990)) are quite difficult. Instead of the microscopic models, mesoscopic coarse-grained models have been proposed and utilized to study the long time relaxation behavior of entangled polymer systems. Due to the nature of the entanglement, however, the mesoscopic models are mainly constructed as phenomenological dynamical models. There are many mesoscopic phenomenological models which reproduce characteristic relaxation behavior of entangled polymers. The Doi-Edwards tube modelDoi and Edwards (1986) is based on the tube picture in which the dynamics of an entangled polymer chain is constrained by a tube like obstacle. The slip-linkHua and Schieber (1998); Masubuchi et al. (2001); Schieber (2003); Doi and Takimoto (2003); Nair and Schieber (2006); Khaliullin and Schieber (2009) and slip-springLikhtman (2005); Uneyama (2011); Chappa et al. (2012); Uneyama and Masubuchi (2012) models employ dynamic links which constrain the motion of polymer chains. The properties of these mesoscopic models depend on the details of the model, and various models and their properties have been studied. Among various mesoscopic coarse-grained models, the multi-chain slip-spring model has some interesting propertiesUneyama and Masubuchi (2012). The multi-chain slip-spring model employs the slip-springs as the extra thermodynamic degrees of freedom, and it has a well-defined (effective) free energy of the system. Both the static and dynamic properties of the model are designed to be statistical mechanically sound, based on the (effective) free energy.

Most of the mesoscopic coarse-grained models employ the entanglement segment as the basic kinetic unit. The coarse-graining levels of these coarse-grained models are similar. Because the numerical costs in simulations depend on the coarse-graining level, highly coarse-grained models with much larger basic kinetic units are preferred to perform simulations for well-entangled systems with very long relaxation times. Kindt and Briels proposed a highly coarse-grained model which is called the responsive particle dynamics (RaPiD) modelKindt and Briels (2007). In the RaPiD model, one polymer chain is simply expressed by a single particle. Compared with other mesoscopic models for entangled polymers, the characteristic length scale of the RaPiD model is large and thus the RaPiD model is numerically efficient. Due to its highly coarse-grained nature, unlike other mesoscopic models, the entanglement effect cannot be directly expressed in the RaPiD model by some mesoscopic objects such as tubes and slip-links. Instead of the tubes and slip-links, Kindt and Briels introduced the transient potential between particles. The transient potential depends on the number of entanglements between two particles, and the number of entanglements is treated as the extra degrees of freedom in the system. In equilibrium, the equilibrium number of entanglements between particles is assumed to be a function of the distance between two particles. Then, by employing an effective free energy for the particle positions and the numbers of entanglements, the dynamics of the system can be modeled by the Langevin equations for the positions and the numbers of entanglements. Kindt and Briels constructed the dynamic equations of the RaPiD model and showed that the coarse-grained simulations with the RaPiD model can successfully reproduce the dynamics of entangled polymers such as the linear viscoelasticity. The RaPiD model has been extended to other systems such as polymer solutionsSantos de Oliveira et al. (2014), associative telechelic polymersSprakel et al. (2009, 2011), and star polymersLiu et al. (2013, 2014); Fitzgerald et al. (2014); Fitzgerald and Briels (2018), and has been shown to be a useful coarse-grained model.

There are some similarities between the multi-chain slip-spring model and the RaPiD model, although their coarse-graining levels are different. Both models employ the extra degrees of freedom to express the entanglement effect, and the effective free energies are utilized to characterize the equilibrium probability distribution. From the similarity between these models, we may consider the RaPiD model as a highly coarse-grained version of the multi-chain slip-spring model, and unify these models. In this work, we propose a RaPiD type highly coarse-grained transient bond model, based on the framework of the multi-chain slip-spring model. We show that we can construct a highly coarse-grained model with the transient bonds and the effective free energy. In our model, transient bonds modulate the dynamics of particles but they do not affect the equilibrium properties. As an ideal case, we consider the system where the equilibrium properties reduce to those of an ideal gas. This ideal version of the transient bond model exhibits interesting dynamical properties even while statically it is an ideal gas. We perform simulations for ideal transient bond systems and study how model parameters affect the dynamic properties such as the diffusion and linear viscoelasticity. We also consider the relation between parameters for the entangled polymer systems and those for the ideal transient bond model.

II Model

II.1 Transient Bond Model

Based on the idea by Kindt and BrielsKindt and Briels (2007), we express one polymer chain by one particle and introduce a transient potential between particles. Unlike the original RaPiD model by Kindt and Briels, we express the transient potential in terms of the discrete number of bonds between particles. We assume that the transient bonds modulate the dynamics of particles but do not affect the equilibrium probability distribution of the particle positions. As we show in what follows, the assumption that a number of transient bonds is discrete is convenient when we consider static, equilibrium properties.

We start from the equilibrium probability distribution without transient bonds. We consider a three dimensional system which consists of $M$ particles and has a volume $V$. We express the position of the $i$-th particle is given as $\bm{R}_{i}$ ($i=1,2,\dots,M$). We assume that the interaction between particles is expressed by two-body interaction potential $v(\bm{r})$. Then the partition function can be simply expressed as

[TABLE]

where $\Lambda$ is the thermal de Broglie wavelength, $\int d\{\bm{R}_{i}\}\dotsb$ represents the integral over all particle positions, $k_{B}$ is the Boltzmann constant, and $T$ is the temperature. The equilibrium probability distribution for the particle positions is simply given as

[TABLE]

To express the dynamic constraint effects, we introduce the transient bonds to the system. We express the number of bonds between the $i$-th and $j$-th particles as $n_{ij}$. (In this model, we assume that multiple bonds can share the same site without any penalties. In other words, we assume transient bonds to be bosons. For convenience, we also assume that bonds at the same site are basically indistinguishable.) Then, the state of the system can be expressed by two sets of variables $\{\bm{R}_{i}\}$ and $\{n_{ij}\}$. We express the transient interaction energy per bond as $u(\bm{r})$ ($\bm{r}$ is the distance vector between particles). Then, the total transient interaction energy becomes

[TABLE]

We should introduce something to control the bond number. In this work, we introduce the effective chemical potential. The idea of the effective chemical potential was first introduced to a slip-link model for entangled polymers by SchieberSchieber (2003), and later applied to the multi-chain slip-spring modelUneyama and Masubuchi (2012). We express the effective chemical potential for a transient bond as $\mu$.

The introduction of the transient bonds affect the equilibrium probability distribution for particle positions, and particles effectively feel an attractive potential. This is the same as the case of the multi-chain slip-spring model, where the slip-springs generate the effective attraction between chainsUneyama and Horio (2011); Chappa et al. (2012); Uneyama and Masubuchi (2012). Such an attractive interaction is an artifact of the model, and should be cancelled so that transient bonds do not affect the equilibrium probability distribution for particles as given by eq (2). For this purpose, we introduce the repulsive compensation potential to the system. We require the following condition for the joint probability distribution of the particle position and the bond number:

[TABLE]

where $P_{\text{eq}}(\{\bm{R}_{i}\})$ should be exactly the same as one given by eq (2), and the summation over $n_{ij}$ is taken for all the possible bond numbers. (Notice that, in general, $P_{\text{eq}}(\{\bm{R}_{i}\})$ can be different from eq (2) if we introduce the bond potentialUneyama and Horio (2011).)

Then, the equilibrium conditional probability distribution should be expressed as follows:

[TABLE]

where $\Xi(\{\bm{R}_{i}\})$ is the grand partition function under given particle positions. (This grand partition is not the (full) grand partition function of the system. The grand partition function in eq (5) should be interpreted as the normalization factor.) The explicit form of the grand partition function becomes

[TABLE]

with the effective fugacity (activity) $\xi\equiv e^{\mu/k_{B}T}$.

The form of the compensation potential is automatically determined by the condition (4). From eqs (2), (5), and (6), the joint equilibrium probability distribution becomes as follows:

[TABLE]

The last term in the exponent of eq (7), which originates from the condition  (4), can be interpreted as the pairwise repulsive compensation potential between particles. The compensation potential cancels the effective attractive interaction by transient bonds, and recovers eq (2) when we take the statistical average over the bond numbers. The effective free energy of the system is defined from eq (7) as

[TABLE]

To simulate dynamical properties, we need dynamic equations. The dynamic equations should satisfy the detailed balance condition in equilibrium. For the dynamics of the particle positions, we employ the overdampled Langevin equation:

[TABLE]

where $\zeta$ is the friction coefficient for a particle, $\bm{\kappa}(t)$ is the velocity gradient tensor, and $\bm{w}(t)$ is the Gaussian white noise. The first and second moments of the noise are given as

[TABLE]

where $\langle\dots\rangle$ represents the statistical average and $\bm{1}$ is the unit tensor. For the dynamics of the transient bonds, we employ a simple birth-death type dynamicsvan Kampen (2007). We assume that the transient bonds are destroyed by a constant rate, and there is no direct correlation between the destruction of different bonds. Then, we have the following transition rate for the decreasing process (the destruction rate):

[TABLE]

where $\tau$ is a constant which represents the average life time of the transient bond. The transition rate for the increasing process (the construction rate) is automatically determined from the detailed balance condition:

[TABLE]

From eqs (7), (11), and (12), we have

[TABLE]

The dynamics of the system can be completely described by eqs (9), (11), and (13). We can perform dynamics simulations by discretizing them and solving simulating time-evolution with some numerical schemes.

Before we proceed to the analysis of the model, we briefly comment on the relation between our model and the original RaPiD model. Our transient bond model may look rather different from the standard RaPiD model. However, if the transient bond number is large, we can approximate the bond number as a continuum variable. Then we can show that our model reduces to the RaPiD type model, as shown in Appendix A. Roughly speaking, the differences between our model and the original RaPiD model are the factor of the bond potential in the effective free energy, and the friction coefficients for particle positions and bond numbers. Thus we believe that our model can be utilized as a simplified version of the RaPiD model.

To calculate the viscoelastic properties, we need an expression of the stress tensor. In many cases, the stress tensor of a mesoscopic models for polymer is assumed to obey by the stress-optical rule. However, in the transient bond model (and also in the RaPiD model), we do not have any intra-chain, conformational degrees of freedom. Thus a naive application of the stress-optical rule does not work. Instead of the stress-optical rule, we simply employ the Kramers form stress tensor:

[TABLE]

In the slip-spring model, the stress tensor consists of two parts; the real stress and the virtual stressRamirez et al. (2007); Uneyama (2011); Uneyama and Masubuchi (2012). The real stress comes from the bond vectors and is consistent with the stress-optical ruleDoi and Edwards (1986); Inoue and Osaki (1996). The virtual stress comes from the slip-springs, and sometimes it is neglected because it is (at least apparently) not consistent with the stress-optical rule. Whether we should include the contribution of the virtual stress or not is not trivialUneyama (2011); Uneyama and Masubuchi (2012). In some cases, the virtual stress term is simply discarded to recover the stress-optical rule. (In the slip-spring model, the contribution of the virtual stress seems not to be so important, at least qualitativelyUneyama (2011); Uneyama and Masubuchi (2012).) In our model, due to the lack of the intra-chain degrees of freedom, we need to include the virtual stress, otherwise the effect of the transient bonds to the viscoelastic properties will be ignored.

Although eq (14) is not (at least apparently) consistent with the stress-optical rule, we should mention that the stress-optical rule is an empirical rule and it does not hold under some situations such as under fast extensional flowsKröger et al. (1997). The simulation data by the Kremer-Grest model showed that the contributions from the non-bonded interaction to the total stress is not simpleRamirez et al. (2007).

The time evolution of the system can be formally expressed by the Fokker-Planck operatorvan Kampen (2007) associated with the Langevin equation (eq (9)) and the reconstruction rates (eqs (11) and (13)). Eq (9) contains the contribution of the external flow, which is absence in equilibrium, and thus we decompose the time evolution operator into the equilibrium and flow parts:

[TABLE]

Here, $\mathcal{L}_{0}$ is the equilibrium time-evolution operator (which is the sum of the Fokker-Planck operator in equilibrium and the transition matrices), and $\mathcal{L}_{1}(t)$ is the time-evolution operator by the external flow:

[TABLE]

If the contribution of the external flow is sufficiently small, the time evolution operator by the flow, $\mathcal{L}_{1}(t)$, can be interpreted as a perturbation. Then we can expand the probability distribution into the perturbation series to consider the linear response:

[TABLE]

By substituting eqs (16) (17) into eq (15), and retaining only the leading order terms, we have

[TABLE]

Here we have utilized $\mathcal{L}_{0}P_{\text{eq}}=0$. Eq (18) can be formally integrated as

[TABLE]

and the ensemble average of the stress tensor at time $t$ is calculated from eq (19):

[TABLE]

where $\langle\dots\rangle_{\text{eq}}$ represents the equilibrium ensemble average, and $\hat{\bm{\sigma}}(t-t^{\prime})\equiv e^{(t-t^{\prime})\mathcal{L}_{0}^{\dagger}}\hat{\bm{\sigma}}$ ($\mathcal{L}_{0}^{\dagger}$ is the adjoint operator for $\mathcal{L}_{0}$) is the time-shifted stress tensor. Eq (20) means that the relaxation modulus tensor is given as the equilibrium auto-correlation function of the stress tensor. For example, from eq (20), the shear relaxation modulus $G(t)$ becomes

[TABLE]

This is consistent with the well-known Green-Kubo formula in the linear response theoryEvans and Morris (2008). (In the slip-spring model, the simple Green-Kubo type form does not hold when we ignore the contribution of the virtual stressRamirez et al. (2007); Uneyama (2011); Uneyama and Masubuchi (2012).)

II.2 Ideal Transient Bond Model

So far, the interaction potential $v(\bm{r})$ and the transient bond potential $u(\bm{r})$ are rather arbitrary. In this work, to investigate dynamical properties of the transient bond model, we set the interaction potential $v(\bm{r})$ to zero and employ a simple harmonic bond potential for $u(\bm{r})$:

[TABLE]

where $a$ is the average bond size. In addition, we limit ourselves to equilibrium systems and set $\bm{\kappa}(t)=0$.

We may call the transient bond model with eqs (22) and (23) as the ideal transient bond model. Eq (22) means that, statically, our system is just an ideal gas from the view point of the static properties in equilibrium. In fact, the partition function of the system is simply calculated to be

[TABLE]

and this gives the equation of state of an ideal gas. (The contributions from the transient bonds and compensation potentials cancel each other.) Therefore, all the static properties of this ideal transient model coincide to those of the ideal gas. For example, the pressure of the ideal transient model in equilibrium is $P=Mk_{B}T/V$. Of course, this does not mean that the dynamical properties of the ideal transient model coincide to the ideal gas. The dynamics of the particles is largely affected by the transient bonds, and thus the dynamical properties such as the viscoelasticity of the ideal transient bond model become qualitatively different from those of the ideal gas without any transient interactions. This would be clear from the Langevin equation (9), where the forces from the bond and compensation potentials do not cancel.

By substituting eq (23) into eq (8), we have the following simple effective free energy:

[TABLE]

Also, the stress tensor (eq (14)) is simplified as

[TABLE]

From eqs (7) and (26), the average stress tensor in equilibrium becomes

[TABLE]

As expected, the contribution of the transient bonds to the stress tensor vanishes in equilibrium, and the stress tensor simply consists of the pressure of the ideal gas.

We can analytically calculate some properties of the ideal transient model. The average number of transient bonds per pair is calculated as

[TABLE]

Thus the transient bond density (the number of transient bonds per unit volume), $\phi$, is given as

[TABLE]

where $\rho=M/V$ is the number density of particles, and we have assumed that the number of particles in the system is sufficiently large ($M\gg 1$). Eq (29) is consistent with the expression for the inter-chain slip-spring density in the multi-chain slip-spring modelUneyama and Masubuchi (2012). Eq (29) means that the transient bond density is proportional to the effective fugacity $\xi$ and the square of the density $\rho$. We can utilize both $\xi$ and $\rho$ to control the transient bond density. (This situation is similar to but much simpler than the multi-chain slip-spring model.) The elastic modulus of the system at the short time scale, $G_{0}$, is roughly proportional to the transient bond density, $G_{0}{\propto\phi}\propto\xi\rho^{2}$.

In the ideal transient bond model, we have only three dimensionless parameters which characterize a target system; the number density of particles $\rho a^{3}$, the effective fugacity $\xi$, and the (relative) average life time of the transient bond $\tau/\tau_{\text{TB}}$ (with $\tau_{\text{TB}}\equiv a^{2}\zeta/k_{B}T$ being the characteristic time scale of the diffusive motion of a particle). For the sake of simplicity, we employ dimensionless units by setting $a=1$, $k_{B}T=1$, and $\zeta$ = 1. Then, the dynamical properties of the system can be fully specified by the set of three parameters, $\rho$, $\xi$, and $\tau$. In what follows we use the dimensionless units and study the effects of these three parameters.

II.3 Parameters for Entangled Polymer Systems

We consider the relation between the degree of polymerization (or the molecular weight) and the dimensionless parameters in the ideal transient bond model, $\rho$, $\xi$, and $\tau$. Because one polymer chain is expressed as a single particle, the characteristic length scale of the ideal transient bond model depends on the degree of polymerization. Due to this property, the relation between an entangled polymer system and the ideal transient bond model is not trivial.

We consider entangled polymer systems where the number of segment density, $\rho_{\text{segment}}$, is constant and only the number of segments (the degree of polymerization), $N$, is changed. We express the size of a segment as $b$. The particle density in the transient bond model can be interpreted as the polymer chain density, and is simply expressed as

[TABLE]

In the ideal transient bond model, we need the dimensionless particle density $\rho a^{3}$ rather than the (dimensioned) particle density itself. Thus we need to specify the average length of the bond $a$. Unfortunately the explicit form of $a$ in terms of polymer parameters is not clear. The transient bonds originate from the entanglement effect, and the entanglement effect becomes relevant only if two chains are overlapped. Then, $a$ should depend on the radius of gyration of a polymer or the average end-to-end distance of a polymer. Fortunately, both of these length scales are proportional to $Nb^{2}$. Therefore it would be reasonable to consider

[TABLE]

In many cases, the unit length scale of a model is taken to be $b$. However, in our model, we should employ $a$ as the unit length scale, because a polymer chain is expressed as a single particle. Therefore, the unit length scale depends on the degree of polymerization rather strongly, $a\propto N^{1/2}$. In other words, the degree of coarse-graining depends on the degree of polymerization rather strongly. This $N$-dependent degree of coarse-graining makes the relation between parameters in the transient bond model and an entangled polymer system rather complicated. Here we consider the $N$-dependence of several physical quantities to establish the relation between transient bond model parameters and the degree of polymerization $N$.

The dimensionless particle density is given as

[TABLE]

Thus the particle density in the transient bond model should be changed if we change the degree of polymerization. The transient bond density can be interpreted as the entanglement density. The (dimensioned) entanglement density is constant if the polymers are sufficiently long and the segment density is constant, $\phi\propto N^{0}$. On the other hand, from eq (29), the dependence of the bond density to the segment number is given as

[TABLE]

To satisfy the condition $\phi\propto N^{0}$ with eq (33), the effective fugacity should depend on $N$. Thus we have the following relation for the effective fugacity.

[TABLE]

The characteristic bond life time $\tau$ is roughly estimated to be the same as the (pure) reptation time $\tau_{\text{rep}}\propto N^{3}$. Also, the (dimensioned) friction coefficient for a particle is related to the friction coefficient for the Rouse chain which moves along the tube, ${\zeta\propto}N\zeta_{\text{segment}}$ ($\zeta_{\text{segment}}$ is the friction coefficient of a segment). Thus we have the following relation for the dimensionless bond life time $\tau/\tau_{\text{TB}}$:

[TABLE]

From the results shown above, all the three parameters in our transient bond model can be approximately related to the degree of polymerization $N$. If we have a reference parameter set in the dimensionless units ($a=1$, $k_{B}T=1$, and $\zeta=1$) as $\rho_{0},\xi_{0}$, and $\tau_{0}$, for the reference degree of polymerization $N_{0}$, we have

[TABLE]

for a system with the degree of polymerization $N$. To convert the dimensionless units in the transient bond model to the standard units for an entangled polymer system, we need the expressions of the dimensionless units in the ideal transient bond model. The unit energy $k_{B}T$ is common for two models, and thus we need only the length and time units, $a$ and $\tau_{\text{TB}}$:

[TABLE]

From (37), we find that the $N$-dependence of the characteristic time unit $\tau_{\text{TB}}$ is the same as one of the Rouse time, $\tau_{R}\propto N^{2}$. Thus we may interpret that the coarse-graining in our model is performed so that the characteristic time scale becomes the Rouse time. This is consistent with eq (35) where we have apparently weak $N$-dependence for $\tau$. Eqs (36) and (37) give only the power law exponents for the degree of polymerization, and they do not tell us about the numerical prefactors. To map the results of the transient bond model to other mesoscopic models for entangled polymers, we should phenomenologically determine the numerical scale conversion factors. Once the scale conversion factors are determined for one system, it is straightforward to obtain the scale conversion factors for other systems with different degrees of polymerization.

II.4 Numerical Scheme

To perform simulations, we discretize time $t$ by the step size $\Delta t$, and set $t_{k}=k\Delta t$. We employ the stochastic Runge-Kutta methodHoneycutt (1992) to integrate the Langevin equation (eq (9)). In the stochastic Runge-Kutta scheme, the update of the position from time $t_{k}$ to $t_{k+1}$ obeys:

[TABLE]

where $\bm{R}_{i,k}\equiv\bm{R}_{i}(t_{k})$, $\bm{F}_{i}$ is the force acting on the $i$-th particle, and $\bm{\theta}_{i,k}$ is the Gaussian white noise. The force is calculated as

[TABLE]

and the noise is a Gaussian noise which is generated to satisfy the following relations:

[TABLE]

We handle each transient bond separately, rather than directly handle the number of the transient bond $n_{ij}$. From eq (11), a transient bond connected to the $i$-th and $j$-th particles will be destroyed by the destruction rate $1/\tau$. The construction rate (eq (13)) is not changed even if we handle existing transient bonds separately. Then, we integrate the construction and destruction rates (eq (11) and $1/\tau$) from $t_{k}$ to $t_{k+1}$, to obtain the destruction probability for an existing bond and the construction probability for a new bond connected to the $i$-th and $j$-th particles:

[TABLE]

Since eq (43) decays rapidly as the distance between particles increases, we assume that the construction of bonds occurs only for pairs of which distance is rather short. Therefore, practically, the construction trials are required only for the pairs whitin a cut-off distance.

The numerical scheme for the dynamics simulation based on the transient bond model is summarized as follows:

Initialize the particle positions and bonds. The particle positions are sampled from the uniform distribution. The bond numbers are sampled from the Poisson distribution, for each pair. 2. 2.

Construct a cell-list for the calculation of the force by eq (40) and the destruction probability by eq (42). 3. 3.

Integrate the Langevin equation for the positions, by the stochastic Runge-Kutta scheme (eqs (38) and (39)). The force is calculated by eq (40) only for the pairs within a cut-off distance. 4. 4.

Destroy the existing bonds by the destruction probability (eq (42)). 5. 5.

Construct new bonds by the construction probability (eq (43)). The construction trials are performed only for the pairs within a cut-off distance. 6. 6.

Go to 2. and iterate the time evolution.

It should be noted here that our numerical scheme shown above is similar to but much simpler than one for the multi-chain slip-spring modelUneyama and Masubuchi (2012). In the slip-spring model, we need to stochastically hop (slide) the slip-springs on polymer chains. We also need to stochastically sample segments around the chain ends for efficient calculations. These are absent in the transient bond model.

During the simulation, the snapshots of the particle positions and the stress tensor of the system are saved. We calculate the average mean-square displacements of particles from the snapshots of particle positions, to study the diffusion behavior. Also, we calculate the shear relaxation moduli by eq (21), to study the viscoelastic behavior. To improve the statistical accuracy, we utilize the Likhtman’s formulaLikhtman (2012) instead of eq (21):

[TABLE]

where $\hat{N}_{\alpha\beta}=\hat{\sigma}_{\alpha\alpha}-\hat{\sigma}_{\beta\beta}$ is the normal stress difference.

III Results

III.1 Effect of Dimensionless Parameters

We perform simulations for different dimensionless parameter sets. To study the effects of individual parameters to the dynamical properties, we employ a parameter set $\rho=1$, $\xi=1$, and $\tau=100$ as a reference parameter set, and change one parameter systematically. We change $\rho$ as $\rho=0.125,0.25,0.5,1$, and $2$ for $\xi=1$ and $\tau=100$, and change $\xi$ as $\xi=0.125,0.25,0.5,1$, and $2$ for $\rho=1$ and $\tau=100$, and change $\tau$ as $\tau=0.1,1,10,100$, and $1000$ for $\rho=1$ and $\xi=1$. The system size is taken to be sufficiently large (typically $V=8^{3}$) and the periodic boundary condition is applied to all the directions. The time step size is set to be $\Delta t=0.01$. The cut-off distance $r_{\text{cut}}$ is set so that $\xi e^{-3r_{\text{cut}}^{2}/2}=10^{-4}$. (Thus the cut-off distance depends on the value of $\xi$.) We perform equilibrium simulations with different random seeds for the same parameter set, and then take averages over time and different samples to improve the statistical accuracy. (We utilize the Mersenne twister methodMatsumoto and Nishimura (1998) to generate random numbers.) We calculate the mean-square displacements and the shear relaxation moduli from the transient bond simulation data.

Fig. 1 shows the relaxation modulus data calculated from the simulation results with different parameter sets. For convenience, the relaxation modulus is normalized by $\rho^{2}\xi$ (which is proportional to the average transient bond density $\phi$, from eq (29)). All the three parameters, $\rho$, $\xi$, and $\tau$ affect the relaxation modulus rather strongly. Especially, the relaxation time increases as these parameters increase. The shapes of the relaxation modulus change as the parameters change, and thus we consider that the transient bond model can reproduce various viscoelastic behavior by tuning the parameters. The longest relaxation time $\tau_{d}$ is estimated from the relaxation modulus at the long time region:

[TABLE]

Fig. 2 shows the dependence of the longest relaxation time to parameters. From Fig. 2(a), we observe that the longest relaxation time strongly depends on $\rho$ and $\xi$. The $\rho$-dependence of the longest relaxation time looks very similar to the $\xi$-dependence. On the other hand, from Fig. 2(b), we observe that the effect of $\tau$ to the longest relaxation time is rather simple. In the large $\tau$ region ($\tau\gtrsim 10$), the longest relaxation time is approximately proportional to the average life time of the bond $\tau$ ($\tau_{d}\propto\tau$, see the dashed line in Fig. 2(b)). This result is physically natural because in our model the stress relaxes when a transient bond is destroyed.

In experiments, the storage and loss moduli measured by an oscillatory shear mode are convenient and widely utilized. We convert the relaxation modulus data into the storage and loss moduli data, by performing the Fourier transform numerically:

[TABLE]

Fig. 3 shows the storage and loss moduli for various parameter sets, calculated from the relaxation modulus data in Fig. 1. As the relaxation modulus, the storage and loss moduli are normalized by $\rho^{2}\xi$. For the cases with large $\rho$ or $\xi$ values ($\rho=2$ and $\xi=2$ in Fig. 3 (a) and (b)), the short time and long time relaxation modes become separated. Other cases exhibit rather broad relaxation mode distributions.

Fig. 4 shows the dependence of the mean-square displacements to parameters $\rho$, $\xi$, and $\tau$. The mean-square displacement decreases as one of the three parameters, $\rho$, $\xi$, and $\tau$, increases. This is consistent with the results for the relaxation modulus. The mean-square displacements exhibit the subdiffusion behavior for relative short time regions. The diffusion coefficient $D$ is estimated from the mean-square displacement at the long time region:

[TABLE]

If there is no transient bonds in the system (the ideal Brownian gas), the diffusion coefficient becomes $D=1$ in the dimensionless unit. Fig. 5 shows the diffusion coefficients estimated from the mean-square displacement data in Fig. 4. For the cases where the parameters $\rho$, $\xi$, or $\tau$ is small, the diffusion coefficient approaches to $D=1$, as expected. As shown in Fig. 5(a), the $\rho$- and $\xi$-dependence of the diffusion coefficient is similar, as the case of the relaxation time. The $\tau$-dependence of the diffusion coefficient is weak. The data in Fig. 5(b) can be fitted to a power-law type relation, $D\propto\tau^{-0.27}$, for the relatively small $\tau$ region, and to constant for the relatively large $\tau$ region. This result seems not to be consistent with the naive expectation from the relaxation time data. We will discuss the diffusion mechanism in the next section.

From these simulation data, we can roughly summarize the behavior of the ideal transient bond model as follows. All the three parameters ($\rho$, $\xi$, and $\tau$) strongly affect the viscoelastic relaxation data, whereas the diffusion data are strongly affected by $\rho$ and $\xi$, and weakly affected by $\tau$.

III.2 Entangled Polymer Systems

We can map the entangled polymers to the ideal transient bond model, when we employ a reference parameter set and use the relations in Sec. II.3. Here we employ $\rho=1$, $\xi=1$, and $\tau=100$ as the reference parameter set for $N=1$. Other simulation conditions are the same as Sec. III.1.

The shear relaxation modulus and the storage and loss moduli for the reference parameter set (in Fig. 1 and Fig. 3) look similar to those of mildly entangled polymers obtained by various mesoscopic simulation models and experiments. In fact, by rescaling the time and stress scales, we can map the storage and loss moduli to those obtained by the Kremer-Grest modelLikhtman et al. (2007) as shown in Fig. 6. The degree of polymerization (number of beads per chain) in the Kremer-Grest model is $N_{\text{KG}}=350$, thus we have the scale conversion factor $N_{\text{TB}}/N_{\text{KG}}=3.5\times 10^{2}$, where $N_{\text{TB}}$ is the degree of polymerization in the transient bond model. The scale conversion factors for the time and modulus (stress) can be determined as the rescaling (or shift) factors used in Fig. 6. The scale conversion factor for the time scales of the transient bond model $t_{\text{TB}}$ and the Kremer-Grest model $t_{\text{KG}}$ is $t_{\text{TB}}/t_{\text{KG}}=8.0\times 10^{3}$, and one for the modulus (stress) scales of the transient bond model $G_{0,\text{TB}}$ and the Kremer-Grest model $G_{0,\text{KG}}$ is $G_{0,\text{TB}}/G_{0,\text{KG}}=3.3\times 10^{-2}$.

We show the shear relaxation modulus data for other values of $N$ ($N=0.25,0.5,2,4$, and $8$) in Fig. 7. The relaxation modulus for large $N$ seems to be sharp. Naively, we consider that this would be because of the lack of the short-time relaxation modes due to the $N$-dependence of the coarse-graining level. Fig. 8 shows the mean-square displacement data for the same values of $N$ as in Fig. 7. The mean-square displacements shown in Fig. 8 correspond to the mean-square displacements of the centers of mass of polymer chains (so-called $g_{3}(t)$ Kremer and Grest (1990)). The mean-square displacement of the center of mass is known to exhibit the crossover behavior Kremer and Grest (1990):

[TABLE]

where $\tau_{e}\propto N^{0}$ and $\tau_{R}\propto N^{2}$ are the entanglement time and the Rouse time, respectivelyDoi and Edwards (1986). The data shown in Fig. 8 are not consistent with eq (49). (For large $N$ cases, we can observe the crossover of the mean-square displacement from constant to the normal diffusion.) To quantitatively analyze the $N$-dependence of the viscoelastic and diffusion behavior, we calculate the longest relaxation time $\tau_{d}$ and the diffusion coefficient $D$. The longest relaxation time $\tau_{d}$ and the diffusion coefficient $D$ of entangled polymer systems are shown in Fig. 9. The longest relaxation time data for relatively large $N$ can be fitted to the power-law, $\tau_{d}\propto N^{3.4}$. This power-law exponent is consistent with the experimental data, the theoretical prediction, and the simulation data by various simulation modelsDoi and Edwards (1986). The diffusion coefficient can be fitted to the power-law, $D\propto N^{-3.3}$, but this power-law exponent is not consistent with experimental data and other simulation modelsLodge (1999).

From these simulation results, we conclude that the viscoelastic behavior of the entangled polymers can be reasonably reproduced if the parameters are determined based on the reference parameter set and the degrees of polymerization. However, the diffusion behavior is not reproduced well. We discuss the effect of the parameters and the mapping of the transient bond model to the entangled polymer systems in detail, in the next section.

IV Discussions

IV.1 Relaxation and Diffusion Behavior of Transient

Bond Model

As shown in Figs. 1-3, the relaxation and diffusion data exhibit similar $\rho$- and $\xi$-dependence. This can be intuitively understood if we consider the average number of transient bonds per particle. From eq (29), the average number of transient bonds per particle is estimated as

[TABLE]

Under the mean-field approximation (which is shown in Appendix B), we expect that the relaxation behavior of a single particle in the system is determined solely by the average number of transient bonds attached to a target particle. Then, the relaxation behavior is approximately determined by the factor $\xi\rho$, and thus the relaxation behavior will be similar if $\xi\rho$ is the same. Conversely, we can almost fully tune the relaxation behavior of the system only by two parameters, $\xi$ and $\tau$, even if we fix $\rho$ to be constant. Such a reduction of the number of free parameters will be useful when we fit the model parameters to a specific target system.

If the value of $\rho\xi$ is sufficiently large, one particle in the system will be strongly constrained by many transient bonds. The relaxation occurs only through the bond reconstruction process, and the characteristic time of the reconstruction is $\tau$. Based on this picture, we expect that the longest relaxation time will approach to $\tau$ for large $\rho\xi$ cases. In Fig. 2(a), we can observe that the longest relaxation time actually approaches to $\tau=100$ for large $\rho$ or large $\xi$ cases. In addition, we expect that the relaxation function will approach to a single exponential form with the relaxation time $\tau$. In other words, the fluctuation of the number of transient bonds attached to one particle will broaden the relaxation mode distribution. Therefore, roughly speaking, we can utilize $\rho$ or $\xi$ to tune the shape of the relaxation mode distribution and utilize $\tau$ to tune the longest relaxation time. These properties are consistent with the estimates under the mean-field approximation in Appendix B.

If the value of $\rho\xi$ is sufficiently small, we expect that the transient bonds cannot form network structures and they can form only dumbbell-like structures (dimers). The longest relaxation time is estimated as one of a dimer (a dumbbell with a harmonic springKröger (2004)), and thus we have $\tau_{d}\approx 1/12$. It should be noted that this relaxation time of a dimer is independent of the bond life time $\tau$. In addition, the contribution of dimers to diffusion coefficient is considered to be small, and thus the diffusion coefficient becomes $D\approx 1$ as we mentioned. This is again independent of the bond life time. Thus the dynamical behavior will be almost independent of the bond life time $\tau$, if $\rho\xi$ is sufficiently small.

The effects of the average life time of the bond $\tau$ on the relaxation time and the diffusion coefficient apparently seem not to be consistent. Unlike the simple $\tau$-dependence of the longest relaxation time, the $\tau$-dependence of the diffusion coefficient seems to be very weak. Moreover, in the large $\tau$ region ($\tau\gtrsim 100$), the diffusion coefficient becomes almost independent of $\tau$. This can be understood as follows. A particle in the system takes two states; the free state in which no bonds are attached to the particle, and the constrained state in which bonds are attached to the particle. In the free state, the diffusion of the particle is not constrained and thus we will observe the free diffusion. In the constrained state, the particle is constrained by the bonds and the average position is almost fixed. The particle at the constrained state can diffuse only via the reconstruction of bonds, thus the diffusion coefficient at the constrained state is inversely proportional to the life time. The average diffusion coefficient is the average of the diffusion coefficients at these two states. For sufficiently large $\tau$, the diffusion coefficient at the constrained state becomes negligibly small and thus we observe the $\tau$-independent diffusion coefficient. (See Appendix B for detailed calculations.) This mechanism is somewhat similar to the dynamic heterogeneity observed in supercooled or glassy systemsYamamoto and Onuki (1998); Sillescu (1999); Uneyama et al. (2015). The transient bond model may be utilized as a model for supercooled or glassy systems which exhibit dynamic heterogeneity.

IV.2 Entangled Polymer Systems

The simulations for the entangled polymer systems by the ideal transient bond model showed reasonable results for the viscoelastic data. The dependence of the relaxation time $\tau_{d}$ to the degree of polymerization $N$ is consistent with the well-known power law, $\tau_{d}\propto N^{3.4}$Doi and Edwards (1986). This result is rather surprising because there is no contour length fluctuation (CLF) in the transient bond model. The pure reptation model gives the power-law exponent $3$, and the exponent $3.4$ is believed to be the apparent exponent due to the correction by the CLF Doi and Edwards (1986). Our result implies that the CLF may not be essential to reproduce the exponent $3.4$. We expect that the fluctuation of the positions of particles may give the correction to the relaxation modulus. Some experimentsLiu et al. ; Matsumiya et al. (2013) and simulationsMasubuchi et al. (2017) report that the CLF mechanism is affected by the constraint-release (CR), and in absence of the CR mechanism, the power-law exponent apparently becomes lower than $3.4$. In the transient bond model, the contribution of the CR type mechanism clearly depends on $N$, and thus the CR might affect the relaxation behavior in a different way from the reptation model. This might be one possible origin of the exponent $3.4$.

The fact that the shear relaxation modulus of the transient bond can be fitted well to that by the Kremer-Grest model is also surprising. The transient bond model is highly coarse-grained, and is designed to reproduce the dynamics at the long time region. Nevertheless, the shear relaxation modulus by the transient bond model agrees well with that by the Kremer-Grest model even at the short time, Rouse relaxation region. Of course, this apparent Rouse like relaxation behavior may be just an artifact. In Fig. 7, we cannot find such Rouse like behavior for systems with large $N$. As we discussed, the relaxation modulus is expected to approach the single exponential form for sufficiently large $N$ because $\rho\xi\propto N$. Then we will only have well-developped plateau at the short time region. Even if the Rouse like behavior is just a model artifact, the reason why we have such a power-law type behavior is not clear. One possible mechanism is that the modulation of the relaxation mode distribution due to the formation of a network-structure. The power-law like viscoelastic behavior at the short time region was experimentally observed in network-forming associative telechelic polymer solutions Uneyama et al. (2012). It is plauseble that the spatial coupling of transient bonds will give similar power-law like relaxation behavior at the short time scale, which (accidentally) has the same power-law exponent as the Rouse model.

In contrast to the viscoelastic behavior, the dependence of the diffusion coefficient $D$ to the degree of polymerization $N$ is not consistent with the well-known power law $D\propto N^{-2.4}$Lodge (1999). We consider that this is due to the lack of the reptation motion in the transient bond model. As we discussed, the diffusion coefficient is determined as the average of the diffusion coefficients at the free and constrained states. The fraction of the free state decreases as $N$ increases. Also, the diffusion coefficients at the constraint states strongly decrease as $N$ increases. Thus we have very strong $N$-dependence of the diffusion coefficient.

We will need to incorporate some mechanisms which reproduce the reptation like diffusion motion to the model, to recover the diffusion behavior which is consistent with experiments and other simulation models. The introduction of another dynamics rule to the model which enhances the diffusion will improve the diffusion behavior. One possible way is to move particles without changing bonds. Such a motion can be realized, for example, if we stochastically exchange the positions of two bonded particles. Other possible ways are to generalize the model and include the conformational degrees of freedomFitzgerald and Briels (2018), and to introduce the reptation type diffusion dynamics by modifying the dynamic equations. Because we have the explicit expression of the effective free energy, such extensions of the model will be rather straightforward.

Although the diffusion behavior of entangled polymers cannot be reproduced well by the transient bond model (at least in the current form), the viscoelastic behavior can be reasonably reproduced. Therefore, we expect that the transient bond model can be utilized to simulate viscoelasticity. Because the transient bond model is highly coarse-grained and is computationally efficient, it will be useful when we are interested only on the viscoelasticity. The relation among our model and other mesoscopic models for entangled polymer systems is interesting. By determining the scale conversion factors among different models, we can connect or compare the simulation data by different modelsMasubuchi and Uneyama (2018). The detailed comparison among our model and some mesoscopic coarse-grained models for entangled polymers is in progress and will be published in near future.

IV.3 Possible Extensions of Model

In this work we performed simulations for the simplest, ideal case of the transient bond model, to investigate the basic properties of the model. We set the interaction potential between particles, $v(\bm{r})$, to be zero, but this is not realistic for polymer melts. The compressiblity of a polymer melt is generally very low, whereas one of the ideal transient bond model is rather high. In the original RaPiD model, Kindt and BrielsKindt and Briels (2007) employed the Gaussian repulsive potential to repel particles. We will also need to employ the Gaussian repulsive potential, to perform more realistic simulations for entangled polymers. However, the entanglement effect exists even at the limit of the zero excluded volume, as far as the chains cannot cross each other. Thus we expect that the ideal transient bond model can capture the characteristic dynamical behavior of entangled polymers even in absence of the Gaussian repulsive potential like the original RaPiD model. It would be worth mentioning here that the multi-chain slip-spring modelUneyama and Masubuchi (2012) can reproduce dynamical properties of entangled polymers even without any interaction potential between segments.

Although in this work we limit ourselves to the ideal systems, the transient bond model can be applied to much complicated systems by using different potentials and dynamics models. Instead of the Langevin equation for the particles, other dynamic equation models can be employed. For example, if we use the dynamic equation of the dissipative particle dynamics (DPD) modelEspañol and Warren (1995); Kinjo and Hyodo (2007), which conserves the momentum, we will be able to simulate the complex flow of viscoelastic materials. Langeloth and coworkers Langeloth et al. (2013) showed that the combination of the multi-chain slip-spring model and the DPD dynamic equation reproduces dynamic properties of entangled polymer melts and solutions reasonably. Due to its high coarse-graining level, the combination of the transient bond model and the DPD dynamic equation will be computationally more efficient than the mutli-chain slip-spring model. The reconstruction dynamics of the transient bonds can be also changed. In some systems, the destruction rate may depend on the bond number or the bond vector. The destruction and construction rates can be changed as long as they satisfy the detailed balance condition. Because our model is based on the effective free energy, such a modification is rather straightforward. The applications of our model to star polymers and telechelic associative polymers will be intersting.

The incorporation of the transient potential Kindt and Briels (2007); Santos de Oliveira et al. (2014); Sprakel et al. (2009, 2011); Liu et al. (2013, 2014); Fitzgerald et al. (2014); Fitzgerald and Briels (2018) is powerful and promising method to model soft matter systems which exhibit complex dynamical behavior. We expect that our approach to unify the RaPiD and slip-spring model can be further generalized, by adding extra degrees of freedom to the system. For example, we can add the average life times of bonds as the extra degrees of freedom, in a similar way to the slip-link modelKhaliullin and Schieber (2009). A recent workUneyama et al. (2015) showed that the diffusion coefficient (or the friction coefficient) should be treated as a fluctuating variable in some systems. In such systems, the diffusion coefficient would be interpreted as extra degrees of freedom to modulate the dynamics, just like the transient bonds in our model. Modeling with extra degrees of freedom to modulate dynamics will be useful for various systems, as an alternative way to the generalized Langevin equation with a memory kernelKawasaki (1973). Our approach would be informative to construct and analyze these dynamical models.

V Conclusion

We constructed the transient bond model, based on the ideas of the RaPiD model and the multi-chain slip-spring model. Our transient bond model has the well-defined effective free energy, and the transient bonds affect only dynamical properties. As the simplest, ideal case, we considered the ideal transient bond model in which the equilibrium properties reduce to those of an ideal gas. The ideal transient bond model can be characterized by three dimensionless parameters, the particle density ${a^{3}}\rho$, the effective fugacity $\xi$, and the average life time of transient bonds $\tau{/\tau_{\text{TB}}}$. The effects of these parameters to the linear viscoelasticity and the mean-square displacement were investigated by simulations. The parameters $\rho$ and $\xi$ have similar effects on the dynamical quantities, because the dynamical behavior is determined by the average number of transient bonds per particle. For entangled polymer systems, we derived the relation between the degree of polymerization and the parameters in the ideal transient bond model. The linear viscoelasticity of the entangled polymer systems can be reasonably reproduced by the ideal transient bond model. However, the diffusion behavior cannot be reproduced by the ideal transient bond model. Thus the transient bond model can be utilized to study entangled polymer systems when we are interested only in the viscoelastic behavior. The transient bond model has a simple structure and it can be tuned for a specific traget system. Also, it can be combined with other mesoscopic models. The application and extension of the transient bond model to much complicated systems will be the future work.

Acknowledgment

The author thanks Prof. Yuichi Masubuchi and Prof. Wim Briels for helpful comments. This work was supported by Grant-in-Aid (KAKENHI) for Scientific Research C 16K05513 and Grant-in-Aid (KAKENHI) for Scientific Research A 17H01152.

Appendix A Approximation for Large Bond Number

In this appendix, we consider the relation between our transient bond model shown in the main text and the (original) RaPiD modelKindt and Briels (2007). In the RaPiD model, the dynamics of the system is described by two Langevin equations. One is the Langevin equation for the position, and is almost the same as that in our model, and another is the Langevin equation for the bond number (in the RaPiD model, the bond number is treated as a continuum variable). The dynamics for the bond number in our model is the birth-death type dynamics, and two models apparently seem to be different.

We consider the case where the number of bonds between the $i$-th and $j$-th particles is large. We consider the dynamics of the bond number $n_{ij}$ under the condition where other variables are fixed. For simplicity, we describe the number of transient bonds as $n$. We express the probability distribution of the bond number $n$ as $P(n,t)$. Then, the dynamic equation for $P(n,t)$ can be expressed as the following master equation:

[TABLE]

where $W(n-1|n)$ and $W(n+1|n)$ are the reconstruction rates (eqs (11) and (13)), and $\bar{n}$ is the average number of the bonds for a fixed bond vector $\bm{R}_{i}-\bm{R}_{j}$,

[TABLE]

We introduce the difference and averaging operators in the bond number space:

[TABLE]

where $f(n)$ is a given function of $n$. With the difference operator $\hat{D}$, eq (51) can be rewritten as follows:

[TABLE]

Here we have defined the flux in the bond number space, $J(n+1/2,t)$:

[TABLE]

Eq (55) has the form of the conservation equation, and if the flux can be related to the difference of the probability distribution, the master equation can be expressed as a Fokker-Planck type equation. The flux can be rewritten as the following form, by utilizing the difference and averaging operators:

[TABLE]

By substituting eq (57) into (55), the master equation can be rewritten as a Fokker-Planck type equationvan Kampen (2007). We assume that $n$ is large, and approximate the difference operator as the differential operator. In addition, we simply ignore the averaging operator. We expand the bond number $n$ around its average value $\bar{n}$, and keep only the leading order terms. Then, the master equation can be approximated as

[TABLE]

Now eq (58) can be interpreted as the Fokker-Planck equation for the bond number $n$. The corresponding Langevin equation for the number of transient bonds between the $i$-th and $j$-th particles, $n_{ij}$, is

[TABLE]

where $\bar{n}_{ij}(\bm{R}_{i}-\bm{R}_{j})\equiv\xi\exp(-u(\bm{R}_{i}-\bm{R}_{j})/k_{B}T)$, and $w_{ij}^{\prime}(t)$ is the Gaussian white noise. The first and second moments of $w_{ij}^{\prime}(t)$ are given as

[TABLE]

The Langevin equation for the particle positions (eq (9)) can be rewritten as follows, in absence of the external flow ($\bm{\kappa}=0$):

[TABLE]

Also, the equilibrium probability distribution (eq (7)) can be approximated by expanding the exponent $n_{ij}$ around the most probable value $\bar{n}_{ij}(\bm{R}_{i}-\bm{R}_{j})$:

[TABLE]

Eqs (61) (59), and (62) have very similar forms to the Langevin equation for the particle positions, one for the numbers of entanglements, and the probability distribution, in the original RaPiD model. (See eqs (3), (6), and (1) in Ref. Kindt and Briels (2007).) In the original RaPiD model, the contribution of the transient bonds to the equilibrium probability distribution is modelled by the harmonic free energy for bond numbers:

[TABLE]

where $\alpha$ is assumed to be constant. We can find a similar factor in the exponent in eq (62). If we set $\alpha$ as $\alpha=1/\bar{n}_{ij}(\bm{R}_{i}-\bm{R}_{j})$ (which is not constant but depends on the bond vector) in eq (63), the equilibrium probability distribution (62) becomes the same as one in the RaPiD model.

From the viewpoint of the dynamics, both our model and the RaPiD model are described by the Langevin equations if the bond number is large. In the RaPiD model, the noise term in the Langevin equation for the particle positions generally depends on the bond vector and the number of entanglements, whereas the noise term in the Langevin equation for the numbers of entanglement is independent of the bond vector. In contrast, the noise term in eq (61) is independent of the bond vector and the bond number. We can interpreted eq (61) as the special case of the RaPiD model, where the friction coefficient is assumed to be constant. On the other hand, the noise term in eq (59) explicitly depends on the bond vector via $\bar{n}_{ij}(\bm{R}_{i}-\bm{R}_{j})$. In general, the Langevin equation for the bond number should be given as the following form, with the effective free energy $\mathcal{F}(\{\bm{R}_{i}\},\{n_{ij}\})$ and the mobility $L_{ij}(\{\bm{R}_{i}\},\{n_{ij}\})$:

[TABLE]

[TABLE]

Comparing eqs (64) and (65) with eq (59), we find $L_{ij}(\{\bm{R}_{i}\},\{n_{ij}\})=\bar{n}_{ij}(\bm{R}_{i}-\bm{R}_{j})/\tau k_{B}T$. Again, we can interpret our model as the special case of the RaPiD model where the friction coefficient for the bond number depends on the equilibrium bond number. Therefore, our model can be interpreted as a variant of the RaPiD model.

Appendix B Mean-Field Approximation For Single Particle Dynamics

To study the dynamical properties of the transient bond model, we consider the mean-field approximation. We consider the statics and dynamics of a single particle in an ideal transient bond model without the external flow. The static property can be obtained by the equilibrium probability distribution for a single particle. Under the mean-field approximation, the information of other particles is smeared out. This situation is almost the same as the relation between the multi-chain slip-spring model and the single-chain slip-spring model as a mean-field model.

We describe the position of a target particle as $\bm{R}$, and the number of transient bonds attached to the target particle as $n$. In the ideal transient bond model, the transient bond is expressed as the harmonic potentials. One end of the $j$-th transient bond is attached to the target particle, and we assume that another end is anchored in space, at $\bm{A}_{j}$. Then, we can approximately express the equilibrium probability distribution as follows:

[TABLE]

where $\tilde{\xi}$ is the effective fugacity, and it does not coincide to the fugacity $\xi$ of the transient bond model in the main text. From eq (66), the mean-field free energy simply becomes the sum of harmonic bond potentials. Thus the Langevin equation can be described as

[TABLE]

where $\bm{w}(t)$ is the Gaussian white noise which satisfies

[TABLE]

The destruction rate of a single transient bond is simply given as

[TABLE]

We assume that the transient bond indices are renumbered after the destruction. The construction rate is determined from the detailed-balance condition. From eqs (66) and (69), we have the construction rate

[TABLE]

where $\bm{A}_{n+1}$ is the position of a new anchoring point.

Some dynamical properties of the transient bond model can be approximately analyzed by the mean-field model shown above. From eqs (69) and (70), the reconstitution process of each transient bond is independent. Thus we have a single exponential type relaxation for the stress:

[TABLE]

The particle is effectively trapped at the average anchoring position unless $n=0$. Thus, for the constrained state ($n\geq 1$), we may rewrite eq (67) as

[TABLE]

with the average anchoring position

[TABLE]

At the short time scale, the particle can freely diffuse but the particle is trapped around the anchoring point, thus at the moderate time scale the motion is strongly constrained. At the long time scale, the particle motion is almost the same as the motion of the average anchoring position ($\bm{R}\approx\bar{\bm{A}}$). Therefore, the mean-square displacement shows the transition from the normal diffusion to subdiffusion, and then to the normal diffusion. The average anchoring position will move roughly by $a/n$ by one reconstruction. If we replace the number of transient bonds by its equilibrium average, we have $a/n\approx a/\tilde{\xi}$. Then, the diffusion coefficient at the constrained state is roughly estimated as

[TABLE]

If the number of bonds is zero, $n=0$, the particle is at the free state and thus it can freely diffuse. The diffusion coefficient at this free state is thus simply estimated as

[TABLE]

The diffusion coefficient is given as the average of the diffusion coefficients at the free and constrained states:

[TABLE]

There are two extreme cases where the diffusion coefficient becomes simple. If the effective fugacity $\tilde{\xi}$ is sufficiently large, the first term in the right hand side of eq (76) becomes negligible and we have $D\approx D_{\text{constrained}}$. If $\tau$ is sufficiently large, the second term in the right hand side of eq (76) becomes negligible and we have $D\approx e^{-\tilde{\xi}}D_{\text{free}}$, and $D$ becomes independent of $\tau$. Although these estimates can not fully explain the simulation results for the transient bond model, we can understand some aspects of the dynamical properties of the transient bond model.