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

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

Takashi Uneyama Center for Computational Science, Graduate School of Engineering, Nagoya University
July 1, 2019
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.

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 particles and has a volume . We express the position of the -th particle is given as (). We assume that the interaction between particles is expressed by two-body interaction potential . Then the partition function can be simply expressed as

 Z=1Λ3MM!∫d{Ri}exp[−1kBT∑i>jv(Ri−Rj)], (1)

where is the thermal de Broglie wavelength, represents the integral over all particle positions, is the Boltzmann constant, and is the temperature. The equilibrium probability distribution for the particle positions is simply given as

 Peq({Ri})=1ZΛ3MM!exp[−1kBT∑i>jv(Ri−Rj)]. (2)

To express the dynamic constraint effects, we introduce the transient bonds to the system. We express the number of bonds between the -th and -th particles as . (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 and . We express the transient interaction energy per bond as ( is the distance vector between particles). Then, the total transient interaction energy becomes

 Utransient({Ri},{nij})=∑i>jniju(Ri−Rj). (3)

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 .

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:

 ∑{nij}Peq({Ri},{nij})=Peq({Ri}), (4)

where should be exactly the same as one given by eq (2), and the summation over is taken for all the possible bond numbers. (Notice that, in general, 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:

 Peq({nij}|{Ri})=1Ξ({Ri})∏i>jnij!exp[μkBT∑i>jnij−1kBT∑i>jniju(Ri−Rj)]. (5)

where 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

 Ξ({Ri})=∏i>j∞∑nij=01nij!exp[μkBTnij−1kBTniju(Ri−Rj)]=exp[ξ∑i>je−u(Ri−Rj)/kBT], (6)

with the effective fugacity (activity) .

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:

 Peq({Ri},{nij})=Peq({nij}|{Ri})Peq({Ri})=1ZΛ3MM![∏i>jξnijnij!]exp[−1kBT∑i>jniju(Ri−Rj)−1kBT∑i>jv(Ri−Rj)−ξ∑i>je−u(Ri−Rj)/kBT]. (7)

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

 F({Ri},{nij})≡∑i>jv(Ri−Rj)+∑i>jniju(Ri−Rj)+kBTξ∑i>je−u(Ri−Rj)/kBT. (8)

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:

 dRi(t)dt=−1ζ∂F({Ri},{nij})∂Ri+κ(t)⋅Ri+√2kBTζwi(t), (9)

where is the friction coefficient for a particle, is the velocity gradient tensor, and is the Gaussian white noise. The first and second moments of the noise are given as

 ⟨wi(t)⟩=0,⟨wi(t)wj(t′)⟩=δijδ(t−t′)1, (10)

where represents the statistical average and 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):

 W(nij−1|nij)=nijτ, (11)

where 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:

 W(nij|nij+1)Peq(nij+1|{Ri})=W(nij+1|nij)Peq(nij|{Ri}). (12)

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

 W(nij+1|nij)=1τξe−u(Ri−Rj)/kBT. (13)

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:

 (14)

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:

 ∂P({Ri},{nij},t)∂t=[L0+L1(t)]P({Ri},{nij},t). (15)

Here, is the equilibrium time-evolution operator (which is the sum of the Fokker-Planck operator in equilibrium and the transition matrices), and is the time-evolution operator by the external flow:

 L1(t)P=−∑i∂∂Ri⋅[κ(t)⋅RiP]. (16)

If the contribution of the external flow is sufficiently small, the time evolution operator by the flow, , can be interpreted as a perturbation. Then we can expand the probability distribution into the perturbation series to consider the linear response:

 P({Ri},{nij},t)=Peq({Ri},{nij})+P1({Ri},{nij},t). (17)

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

 ∂P1({Ri},{nij},t)∂t≈L1(t)Peq({Ri},{nij})+L0P1({Ri},{nij},t)=VkBT^σ:κ(t)Peq({Ri},{nij})+L0P1({Ri},{nij},t). (18)

Here we have utilized . Eq (18) can be formally integrated as

 P1({Ri},{nij},t)=VkBT∫t−∞dt′e(t−t′)L0^σ:κ(t′)Peq({Ri},{nij}), (19)

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

 σ(t)=∫d{Ri}∑{nij}^σ[Peq({Ri},{nij})+P1({Ri},{nij},t)]=⟨^σ⟩eq+VkBT∫t−∞dt′⟨^σ(t−t′)^σ⟩eq:κ(t′), (20)

where represents the equilibrium ensemble average, and ( is the adjoint operator for ) 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 becomes

 G(t)=VkBT⟨^σxy(t)^σxy⟩eq. (21)

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 and the transient bond potential are rather arbitrary. In this work, to investigate dynamical properties of the transient bond model, we set the interaction potential to zero and employ a simple harmonic bond potential for :

 v(r) =0, (22) u(r) =3kBT2a2r2, (23)

where is the average bond size. In addition, we limit ourselves to equilibrium systems and set .

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

 Z=VMΛ3MM!, (24)

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

 F({Ri},{nij})kBT=32a2∑i>jnij(Ri−Rj)2+ξ∑i>je−3(Ri−Rj)2/2a2. (25)

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

 ^σ=3kBT2a2V∑i>j[nij−ξe−3(Ri−Rj)2/2a2](Ri−Rj)(Ri−Rj)−MkBTV1. (26)

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

 ⟨^σ⟩eq=∫d{Ri}∑{nij}^σPeq({nij}|{Ri})Peq({Ri})=−MkBTV1. (27)

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

 ⟨nij⟩eq=∫dRidRj∑nijnijPeq(Ri,Rj,nij)=∫dRidRj1V2∞∑nij=0nijnij!ξnijexp[−3nij(Ri−Rj)22a2−ξe−3(Ri−Rj)2/2a2]=ξV(2πa23)3/2. (28)

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

 ϕ=M(M−1)2V⟨nij⟩eq≈ξρ22(2πa23)3/2, (29)

where is the number density of particles, and we have assumed that the number of particles in the system is sufficiently large (). 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 and the square of the density . We can utilize both and 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, , is roughly proportional to the transient bond density, .

In the ideal transient bond model, we have only three dimensionless parameters which characterize a target system; the number density of particles , the effective fugacity , and the (relative) average life time of the transient bond (with being the characteristic time scale of the diffusive motion of a particle). For the sake of simplicity, we employ dimensionless units by setting , , and = 1. Then, the dynamical properties of the system can be fully specified by the set of three parameters, , , and . 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, , , and . 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, , is constant and only the number of segments (the degree of polymerization), , is changed. We express the size of a segment as . The particle density in the transient bond model can be interpreted as the polymer chain density, and is simply expressed as

 ρ=ρsegmentN. (30)

In the ideal transient bond model, we need the dimensionless particle density rather than the (dimensioned) particle density itself. Thus we need to specify the average length of the bond . Unfortunately the explicit form of 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, 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 . Therefore it would be reasonable to consider

 a2∝Nb2. (31)

In many cases, the unit length scale of a model is taken to be . However, in our model, we should employ 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, . In other words, the degree of coarse-graining depends on the degree of polymerization rather strongly. This -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 -dependence of several physical quantities to establish the relation between transient bond model parameters and the degree of polymerization .

The dimensionless particle density is given as

 ρa3∝ρsegmentNN3/2b3∝N1/2. (32)

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, . On the other hand, from eq (29), the dependence of the bond density to the segment number is given as

 ϕ=ξρ22(2πa23)3/2∝ξρ2segmentN2N3/2b3∝ξN−1/2. (33)

To satisfy the condition with eq (33), the effective fugacity should depend on . Thus we have the following relation for the effective fugacity.

 ξ∝N1/2. (34)

The characteristic bond life time is roughly estimated to be the same as the (pure) reptation time . Also, the (dimensioned) friction coefficient for a particle is related to the friction coefficient for the Rouse chain which moves along the tube, ( is the friction coefficient of a segment). Thus we have the following relation for the dimensionless bond life time :

 ττTB=kBTτa2ζ∝kBTτrepNb2×Nζsegment∝N. (35)

From the results shown above, all the three parameters in our transient bond model can be approximately related to the degree of polymerization . If we have a reference parameter set in the dimensionless units (, , and ) as , and , for the reference degree of polymerization , we have

 ρ=(NN0)1/2ρ0,ξ=(NN0)1/2ξ0,τ=(NN0)1τ0, (36)

for a system with the degree of polymerization . 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 is common for two models, and thus we need only the length and time units, and :

 a∝N1/2b∝N1/2,τTB=a2ζkBT∝N2b2ζsegmentkBT∝N2. (37)

From (37), we find that the -dependence of the characteristic time unit is the same as one of the Rouse time, . 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 -dependence for . 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 by the step size , and set . 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 to obeys:

 R∗i,k =Ri,k+ΔtFi({Ri,k})+√2Δtθi,k, (38) Ri,k+1 =Ri,k+Δt2[Fi({Ri,k})+Fi({R∗i,k})]+√2Δtθi,k, (39)

where , is the force acting on the -th particle, and is the Gaussian white noise. The force is calculated as

 Fi({Ri})≡−3∑j[nij−ξe−3(Ri−Rj)2/2](Ri−Rj), (40)

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

 ⟨θi,k⟩=0,⟨θi,kθj,l⟩=δijδkl1. (41)

We handle each transient bond separately, rather than directly handle the number of the transient bond . From eq (11), a transient bond connected to the -th and -th particles will be destroyed by the destruction rate . 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 ) from to , to obtain the destruction probability for an existing bond and the construction probability for a new bond connected to the -th and -th particles:

 Ψ− =1−exp(−Δt/τ), (42) Ψ+,ij =1−exp[−ξe−3(Ri−Rj)2/2Δt/τ]. (43)

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:

1. 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. Construct a cell-list for the calculation of the force by eq (40) and the destruction probability by eq (42).

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. Destroy the existing bonds by the destruction probability (eq (42)).

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. 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):

 G(t)=V5kBT[⟨^σxy(t)^σxy⟩eq+⟨^σyz(t)^σyz⟩eq+⟨^σxz(t)^σzx⟩eq]+V30kBT[⟨^Nxy(t)^Nxy⟩eq+⟨^Nyz(t)^Nyz⟩eq+⟨^Nzx(t)^Nzx⟩eq], (44)

where 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 , , and as a reference parameter set, and change one parameter systematically. We change as , and for and , and change as , and for and , and change as , and for and . The system size is taken to be sufficiently large (typically ) and the periodic boundary condition is applied to all the directions. The time step size is set to be . The cut-off distance is set so that . (Thus the cut-off distance depends on the value of .) 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 (which is proportional to the average transient bond density , from eq (29)). All the three parameters, , , and 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 is estimated from the relaxation modulus at the long time region:

 lnG(t)≈(const)−τdt(t≳τ). (45)

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 and . The -dependence of the longest relaxation time looks very similar to the -dependence. On the other hand, from Fig. 2(b), we observe that the effect of to the longest relaxation time is rather simple. In the large region (), the longest relaxation time is approximately proportional to the average life time of the bond (, 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:

 G′(ω) =ω∫∞0dtG(t)sin(ωt), (46) G′′(ω) =ω∫∞0dtG(t)cos(ωt). (47)

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 . For the cases with large or values ( and 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 , , and . The mean-square displacement decreases as one of the three parameters, , , and , 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 is estimated from the mean-square displacement at the long time region:

 ln⟨[Ri(t)−Ri(0)]2⟩≈ln(6D)+lnt(t≳τ). (48)

If there is no transient bonds in the system (the ideal Brownian gas), the diffusion coefficient becomes 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 , , or is small, the diffusion coefficient approaches to , as expected. As shown in Fig. 5(a), the - and -dependence of the diffusion coefficient is similar, as the case of the relaxation time. The -dependence of the diffusion coefficient is weak. The data in Fig. 5(b) can be fitted to a power-law type relation, , for the relatively small region, and to constant for the relatively large 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 (, , and ) strongly affect the viscoelastic relaxation data, whereas the diffusion data are strongly affected by and , and weakly affected by .

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 , , and as the reference parameter set for . 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 , thus we have the scale conversion factor , where 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 and the Kremer-Grest model is , and one for the modulus (stress) scales of the transient bond model and the Kremer-Grest model is .

We show the shear relaxation modulus data for other values of (, and ) in Fig. 7. The relaxation modulus for large seems to be sharp. Naively, we consider that this would be because of the lack of the short-time relaxation modes due to the -dependence of the coarse-graining level. Fig. 8 shows the mean-square displacement data for the same values of 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 Kremer and Grest (1990)). The mean-square displacement of the center of mass is known to exhibit the crossover behavior Kremer and Grest (1990):

 ⟨[Ri(t)−Ri(0)]2⟩∝⎧⎪⎨⎪⎩t1(t≲τe),t1/2(τe≲t≲τR),t1(τR≲t), (49)

where and 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 cases, we can observe the crossover of the mean-square displacement from constant to the normal diffusion.) To quantitatively analyze the -dependence of the viscoelastic and diffusion behavior, we calculate the longest relaxation time and the diffusion coefficient . The longest relaxation time and the diffusion coefficient of entangled polymer systems are shown in Fig. 9. The longest relaxation time data for relatively large can be fitted to the power-law, . 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, , 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 - and -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

 ϕρ=ξρ2(2πa33)3/2∝ξρ. (50)

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 , and thus the relaxation behavior will be similar if is the same. Conversely, we can almost fully tune the relaxation behavior of the system only by two parameters, and , even if we fix 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 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 . Based on this picture, we expect that the longest relaxation time will approach to for large cases. In Fig. 2(a), we can observe that the longest relaxation time actually approaches to for large or large cases. In addition, we expect that the relaxation function will approach to a single exponential form with the relaxation time . 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 or to tune the shape of the relaxation mode distribution and utilize 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 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 . It should be noted that this relaxation time of a dimer is independent of the bond life time . In addition, the contribution of dimers to diffusion coefficient is considered to be small, and thus the diffusion coefficient becomes 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 , if is sufficiently small.

The effects of the average life time of the bond on the relaxation time and the diffusion coefficient apparently seem not to be consistent. Unlike the simple -dependence of the longest relaxation time, the -dependence of the diffusion coefficient seems to be very weak. Moreover, in the large region (), the diffusion coefficient becomes almost independent of . 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 , the diffusion coefficient at the constrained state becomes negligibly small and thus we observe the -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 to the degree of polymerization is consistent with the well-known power law, 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 , and the exponent 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 . 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 . In the transient bond model, the contribution of the CR type mechanism clearly depends on , 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 .

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 . As we discussed, the relaxation modulus is expected to approach the single exponential form for sufficiently large because . 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 to the degree of polymerization is not consistent with the well-known power law 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 increases. Also, the diffusion coefficients at the constraint states strongly decrease as increases. Thus we have very strong -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, , 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 , the effective fugacity , and the average life time of transient bonds . The effects of these parameters to the linear viscoelasticity and the mean-square displacement were investigated by simulations. The parameters and 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 -th and -th particles is large. We consider the dynamics of the bond number under the condition where other variables are fixed. For simplicity, we describe the number of transient bonds as . We express the probability distribution of the bond number as . Then, the dynamic equation for can be expressed as the following master equation:

 ∂P(n,t)∂t=W(n|n+1)P(n+1,t)+W(n|n−1)P(n−1,t)−[W(n+1|n)+W(n−1|n)]P(n,t)=1τ[(n+1)P(n+1,t)+¯nP(n−1,t)−(¯n+n)P(n,t)], (51)

where and are the reconstruction rates (eqs (11) and (13)), and is the average number of the bonds for a fixed bond vector ,

 ¯n≡ξe−u(Ri−Rj)/kBT. (52)

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

 ^Df(n) ≡f(n+1/2)−f(n−1/2), (53) ^Mf(n) ≡f(n+1/2)+f(n+1/2)2, (54)

where is a given function of . With the difference operator , eq (51) can be rewritten as follows:

 ∂P(n,t)∂t=−1τ[[¯nP(n,t)−(n+1)P(n+1,t)]−[¯nP(n−1,t)−nP(n,t)]]=−^DJ(n,t). (55)

Here we have defined the flux in the bond number space, :

 J(n+1/2,t)≡1τ[¯nP(n,t)−(n+1)P(n+1,t)]. (56)

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:

 J(n+1/2,t)=1τ[¯nP(n,t)−(n+1)P(n+1,t)]=1τ[¯n(^M−^D/2)P(n+1/2,t)−(n+1)(^M+^D/2)P(n+1/2,t)]=−1τ[(n+1)+¯n2[2(n+1)−¯n(n+1)+¯n^MP(n+1/2,t)+^DP(n+1/2,t)]]. (57)

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

 (58)

Now eq (58) can be interpreted as the Fokker-Planck equation for the bond number . The corresponding Langevin equation for the number of transient bonds between the -th and -th particles, , is

 dnij(t)dt=−nij−¯nij(Ri−Rj)τ+√2¯nij(Ri−Rj)τw′ij(t), (59)

where , and is the Gaussian white noise. The first and second moments of are given as

 ⟨w′ij(t)⟩=0,⟨w′ij(t)w′kl(t′)⟩=δij,klδ(t−t′). (60)

The Langevin equation for the particle positions (eq (9)) can be rewritten as follows, in absence of the external flow ():

 dRi(t)dt=−1ζ∑j[∂v(Ri−Rj)∂Ri+[nij−¯nij(Ri−Rj)]∂u(Ri−Rj)∂Rj]+√2kBTζwi(t). (61)

Also, the equilibrium probability distribution (eq (7)) can be approximated by expanding the exponent around the most probable value :

 Peq({Ri},{nij})≈1ZΛ3MM![∏i>j1√2π¯nij(Ri−Rj)]×exp[−1kBT∑i>jv(Ri−Rj)−∑i>j[nij−¯nij(Ri−Rj)]22¯nij(Ri−Rj)]. (62)

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:

 Fbond({Ri},{nij})=∑i>jkBTα2[nij−¯nij(Ri−Rj)]2 (63)

where is assumed to be constant. We can find a similar factor in the exponent in eq (62). If we set as (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 . In general, the Langevin equation for the bond number should be given as the following form, with the effective free energy and the mobility :

 dnij(t)dt=−Lij({Ri},{nij})∂F({Ri},{nij})∂nij+kBT∂Lij({Ri},{nij})∂nij+√2kBTL(n)ij({Ri},{nij})w′ij(t), (64)