text:
# Dynamical Properties of the Slithering Snake Algorithm: A numerical test of the activated reptation hypothesis

## I Introduction

#### Background.

The molecular-level description of the dynamics of non-dilute polymer liquids remains a challenging problem of great fundamental interest as well as practical importance LRP90 . In the last 30 years, attention has largely focused on the applicability of the reptation concept introduced by de Gennes and Edwards DG79 ; DE86 . This approach postulates that individual chains, constrained by their neighbors, move primarily along their own contours in a snakelike fashion. Reptation-based models have been very successful in reconciling a wide range of experimental observations LRP90  and are broadly in agreement with more recent computational results as reviewed in KP97 ; KG95 . All evidence amounts to the existence of an additional and chain length independent length scale. This length scale is called the “tube diameter” $`d_\text{e}N_\text{e}^\nu `$ with $`N_\text{e}(\varphi ,l_\text{p})`$ the associated “entanglement length”, $`\varphi `$ the molecular density, $`l_p`$ the persistence length and $`\nu `$ the Flory exponent of the polymer chain. The density and persistence length dependence of $`d_\text{e}`$ have been probed by various numerical methods BB00 ; PBHK91 ; PBHK91\_JCP ; MWB00 ; FM01 , indicating that it is proportional to the excluded volume blob size $`\xi (\varphi ,l_\text{p})`$. While, strictly speaking, not having the (more prestigious) status of a ‘scientific theory’, the reptation concept is widely (though not generally) seen as a physically appealing starting point any further modeling approach has to compete with.

#### First motivation.

The simulation of entangled polymer melts is a challenging task and it is a matter of debate if the (suspected) reptating motion of polymer chains has already been reproduced numerically in a sufficiently large simulation box PKG00 . Without any doubt, it will be near to impossible to simulate within the near future by means of locally realistic dynamics industrially relevant problems, such as the rheology of polymer melts with fillers. It is, therefore, tempting to accept the snakelike motion as a basic ingredient in a coarse-grained algorithm and to verify if a meaningful dynamic interpretation of the well-known “slithering snake algorithm” (SSA) Kron ; Wall ; Sokal  can be given. This algorithm is sketched in Fig. 1. It has been extensively used for three decades in simulations of static properties because of its efficiency in exploring phase space. Clearly, it would be extremely useful if a (simple and transparent) mapping of this algorithm on its corresponding local dynamical scheme could be achieved.

#### Self-consistency of reptation.

While most of the current computational work appears to focus on the demonstration of the curvilinear chain motion KG95 ; MWB00 ; FM01 ; KBM01 ; PKG00 , here we take this key reptation hypothesis for granted. Therefore, we are not interested in the Rouse-like dynamics at short times and for chain length smaller than $`N_\text{e}`$. For the SSA we expect $`N_\text{e}1`$ so that it should be possible to explore the dynamics of the so-called ‘primitive path’ DE86 . What we do question, however, is an additional (and arguably less profound) mean-field assumption made in the original versions of reptation which supposes the uncorrelated motion of neighboring “snakes”.111We will often use the term “snake” as a short form for “chains moving according to the slithering snake algorithm”. In the following we refer to this assumption as the “uncorrelated reptation hypothesis” (URH). One of its immediate consequences is, for instance, the linear time dependence of the mean-square displacements (MSD) of the center of mass or of the curvilinear motion of an inner chain monomer for the snake dynamics. The verification of the influence of the interchain correlations on reptating snakes is the second motivation of this study.

#### Internal modes, tube fluctuations.

In recent years a number of theoretical models and numerical schemes DE86 ; Doi83 ; Rubinstein87 ; DesCloizeaux  have been developed in order to explain the observed molecular weight dependence of the zero-shear viscosity $`\eta N^{3.4\pm 0.2}`$ for $`NN_e`$ LRP90 ; DE86 . The main idea of these theories is to account for the internal modes of polymer dynamics inside the reptation tube (“tube-length fluctuations”) or to introduce a finite lifetime of the tube (“double reptation” DesCloizeaux ) without addressing, however, the possibility of correlated chain motion. While these are undoubtedly valiant efforts, they are of no concern here, as we are, to stress it again, interested in the physics at much larger molecular weights and correspondingly longer time scales.

#### Coupling between dynamics and density fluctuations.

The possibility of a coupling between primitive path dynamics and the fluctuating excluded volume interactions of neighboring chains was first investigated by Deutsch Deutsch . He suggested that the reptation of a chain tail from its original tube to a new environment implies a finite free energy penalty $`F(s)`$ proportional to the curvilinear length of the tail $`s`$. Moreover, it was assumed that this additional activation energy of the snake can relax as soon as the spatial size of the escaped portion becomes of the order of the mean distance between chain ends $`d_{\text{end}}N^{1/3}`$. Taking into account the Gaussian chain statistics in a dense polymer solution, Deutsch obtained an exponential increase $`\tau \eta N^3\mathrm{exp}(N^{2/3})`$ for the relaxation times and viscosities in concentrated polymer solutions with extremely large molecular masses. More recently, these ideas have been elaborated further by Semenov and Rubinstein Semenov91 ; Semenov97 ; SR98  and implemented in a slightly different and thermodynamically fully consistent way. They predict that the URH breaks down for $`NN_\text{e}^3`$ SR98 . Since we expect $`N_\text{e}`$ to be of order 1 for the SSA, it should be possible to test this “activated reptation” hypothesis (ARH) in our simulations.

#### Brief summary of our work.

We have harnessed a SSA of the bond-fluctuation lattice model and compared its dynamics for consistency with the URH of the classical reptation model. We find that the dynamics rapidly slows down with chain overlap and becomes strongly subdiffusive in a finite time interval $`\tau _{}t\tau _+`$ whose bounds are determined by the blob size $`\xi `$ and the mean distance between chain ends $`d_{\text{end}}`$. This slowing down is at variance with the URH and indicates strong correlations between neighboring chains. To corroborate our results for the SSA dynamics in an annealed molecular field, referred to as “annealed SSA” (aSSA) hereafter, we have considered a different version of the SSA where all chains but one have been quenched. We will refer to this version as the “quenched SSA” (qSSA). The comparison of both dynamical schemes demonstrates the mutual influence of the correlation hole retaining the snakes by means of an effective force and the fluctuating chain ends which allow the snakes to sneak out of the local cage of immobile monomer sites.

#### Outline.

This work is presented as follows. Our numerical scheme is briefly recapitulated in the next section. The parameter range and some technical points concerning the configurations and measurements are presented in Sec. III. The relevant analytical predictions according to the URH and ARH are summarized in Sec. IV. The numerical results of this paper are presented in Sec. V where we characterize and compare aSSA and qSSA dynamics. The implications and possible extensions of our work are discussed in the final Sec. VI. There, we speculate on the importance of local motion to obtain a fully consistent mapping of a slithering snake scheme on a locally realistic model.

## II Simulation model

#### Bond-fluctuation model.

In this investigation we use the extensively studied ‘bond-fluctuation model’ (BFM) PBHK91 ; DB91 ; Marcus , a lattice Monte Carlo scheme introduced by Carmesin and Kremer CK88 . Figure 1 depicts a two dimensional projection of both the original BFM with local monomer dynamics (upper chain) and the BFM with global snakelike moves (lower chain). These versions have been discussed elsewhere in detail DB91 ; Marcus ; WBB95  and we concentrate here on some salient features.

In the BFM each monomer occupies the eight corners of the unit cube on a simple cubic lattice. The bonds between adjacent monomers can vary in length and direction, subject only to excluded volume constraints and entanglement restrictions. Basically, due to the larger number of degrees of freedom — $`87`$ angles between bonds are possible — the BFM does not suffer from the same ergodicity problems that classical Monte Carlo schemes of self-avoiding random chains on simple cubic lattices are prone to Wall ; Sokal . (We come back to this point in Sec. V.5.) In the original BFM version only local jumps to nearest neighbor sites are allowed. The excluded volume then automatically forbids any intersection of bonds.

#### Slithering snake schemes: aSSA and qSSA.

If snakelike moves are considered, the uncrossability of bonds has to be checked explicitly to avoid configurations which cannot be attained or unraveled by the original local scheme WP96 . The snake move along the existing backbone of a chain is equivalent to removing an end monomer and connecting it to the other end of the chain, leaving the lattice positions of the middle monomers unchanged. Therefore, density fluctuations occur only at chain ends if no local moves are allowed.

The main focus of this study is the dynamics of the “annealed SSA” (aSSA) where all chains are allowed to move: A chain and one of its ends are selected at random and a global move in a randomly chosen bond direction is attempted. We compare this with the “quenched SSA” (qSSA) of systems where all chains but one are quenched, thus, disallowing temporal density fluctuations (at chain ends or elsewhere). We sample the dynamics of the remaining free snake and average over all chains of one quenched equilibrium configuration. Unless stated otherwise, in the following we always refer to aSSA properties. Data from aSSA (qSSA) are generally indicated by open (full) symbols.

#### Definition of a Monte Carlo step.

Following Sokal Sokal  we have chosen to set the unit of the Monte Carlo time step (MCS) such that on average each monomer is attempted once per cycle, irrespective of whether this attempt occurs through local or snake moves. If only the latter are allowed, every chain is on average visited once per MCS. Our definition of the time scale differs from the one used traditionally in the description of slithering snakes, which focuses on the computational time required to equilibrate samples. Indeed every snake move requires (essentially) the same CPU time irrespective of $`N`$, since all monomers are moved simultaneously. Hence, the computational relaxation time is proportional to $`\tau _\text{t}/N`$ where $`\tau _\text{t}`$ stands for the terminal relaxation time in our units.

## III Parameters, configurations, measurements

This section presents various technical aspects of the simulation and introduces the important static and dynamical quantities for the subsequent analysis. When a quantity is defined, we often refer to the figures where its dependence on $`N`$, $`\varphi `$ or time is shown, without discussing the figure in detail, however. This discussion will be done in Sec. V.

#### Parameters.

In this paper we study athermal systems (no additional potential parameters) containing monodisperse snakes (i.e., disallowing any local move) at volume fractions $`\varphi `$ ranging from dilute solutions to dense melts ($`\varphi 0.5`$). The volume fraction is defined by $`\varphi =8NN_\text{p}/L^3`$, where $`N`$, $`N_\text{p}`$ and $`L`$ are the chain length, the number of polymers in the system and the linear dimension of the simulation box, respectively. In the present work, we primarily investigate the scaling of dynamical properties with mass $`N`$ (see, e.g., Figs. 2,69,17,18). Chain lengths ranging from $`N=16`$ up to $`N=1024`$ have been simulated. This is sufficient to put theoretical predictions to the test. We use a periodic cubic lattice of linear size $`L=128`$. At $`\varphi =0.5`$ this corresponds to $`131072`$ monomers per simulation box. Note that all length scales are given in units of the lattice spacing $`a`$ (i.e., $`a=1`$ in the following).

In fact, we report here the first results from a larger and still on-going study, where in addition to $`N`$ and $`\varphi `$ we have also varied the ratio $`\omega `$ of local to snake moves over several orders of magnitude ranging continuously from purely local ($`\omega =\mathrm{}`$) to pure snake dynamics ($`\omega =0`$). For the sake of comparison, the results of the local dynamics are included in some figures (Figs. 3,18). At the end of this paper, we briefly comment on the influence that a finite value for $`\omega `$ has on the dynamics (Fig. 18).

#### Measurement of static properties.

Obviously, the static properties of an ergodic system do not depend on the dynamics of the algorithm. This has been carefully checked for the annealed SSA. (Numerically, it also appears to be true for the qSSA, even though ergodicity is less obvious in this case, since the environment of reference chains is frozen. Note however that all initial qSSA chains are taken from an equilibrated annealed ensemble.) This allows us to use the precise characterization obtained in previous studies, e.g., for the mean bond length $`l`$, the mean end-to-end distance $`R_\text{e}`$ or the mean radius of gyration $`R_\text{g}`$. The first quantity gives the mean curvilinear distance along the chain associated with every successful snake move, the second the typical displacement $`R_\text{e}/N`$ of the center of mass per snake move. We recall that a polymer chain has a fractal dimension $`1/\nu `$ with $`\nu `$ being the Flory exponent: $`R_\text{e}R_\text{g}N^\nu `$. The Flory exponent for dilute chains in a good solvent is $`\nu =\nu _03/5`$. Above the overlap density $`\varphi ^{}N/R_\text{g}^3N^{13\nu _0}`$, the chains become Gaussian chains of blobs DG79  with $`\nu =1/2`$. An important length scale for the subsequent discussion is the mean distance between chain ends, $`d_{\text{end}}=(4N/\varphi )^{1/3}`$, which we will often refer to as “mean end distance” in the following. Frequently, we use the “effective bond length” DE86  $`b(\varphi )lim_N\mathrm{}R_\text{e}/N^{\nu =1/2}`$ which is weakly density dependent. This definition is chain length independent only for systems above $`\varphi ^{}`$. In the dilute limit, we have $`R_\text{e}/N^{\nu _0}`$, which yields $`b_03`$ for $`N1`$. Some of these properties are summarized in Tables 1 and 2 for $`\varphi =0.5`$ and $`\varphi =0.125`$, respectively. Note that $`R_\text{g}d_{\text{end}}=2N^{1/3}`$ at $`\varphi =0.5`$ for $`N16`$, while for $`\varphi =0.125`$ chains with $`N<64`$ do not overlap, $`R_\text{g}<d_{\text{end}}`$.

#### Various definitions of dynamic properties measured.

As dynamical properties we have sampled various spatial and curvilinear mean-square displacements (MSD) $`g_i(t)`$, as well as associated times, diffusion coefficients and histograms. Following previous work PBHK91  we define the spatial MSD $`g_1(t)`$ for the middle monomer (with monomer index $`n=N/2`$), $`g_2(t)`$ for the motion of the middle monomer in the frame of the center of mass, $`g_3(t)`$ for the center of mass motion, $`g_4(t)`$ for the end monomer motion (averaged over $`n=1`$ or $`n=N`$) and $`g_5(t)`$ for the motion of the end monomers in the center of mass frame. Examples are given in Figs. 45. The MSD related to the curvilinear motion along the backbone is called $`g_6(t)=s(t)^2`$ where $`s`$ is the curvilinear distance traveled by the middle monomer $`n=N/2`$ in time $`t`$. It is presented in Figs. 12, 13 and 14. As we will see below, the scaling of $`g_6(t)`$ and the associated histograms (see Figs. 15,16,17) are crucial for the understanding of the slithering snake dynamics.

The curvilinear and spatial diffusion coefficients, presented in Tables 1 and 2 and in Figs. 6, 7 and 18, are defined by $`D_\text{c}lim_t\mathrm{}g_6(t)/2t`$ and $`D_\text{s}lim_t\mathrm{}g_3(t)/6t`$. From the diffusion coefficients we introduce the time scales $`\tau _\text{s}R_{\text{g}}^{}{}_{}{}^{2}/D_\text{s}`$ and $`\tau _\text{c}\left(Nl\right)^2/(\pi ^2D_\text{c})`$. The prefactor in the last definition has been introduced for convenience (see the discussion of Fig. 15 in Sec. V.5). Additional times can be introduced to characterize the different MSD’s. Following again Ref. PBHK91  the time scales $`\tau _1`$, $`\tau _2`$, $`\tau _3`$ and $`\tau _4`$ are defined by $`g_1(t\tau _1)=R_\text{g}^2`$, $`g_2(t\tau _2)=2R_\text{g}^2/3`$, $`g_2(t\tau _3)=g_3(t=\tau _3)`$ and $`g_5(t\tau _4)=R_\text{g}^2`$. In addition to this we shall use $`g_1(t\tau _{\text{end}})=d_{\text{end}}^2`$ and $`g_1(t\tau _d)=d^2`$ where $`d`$ is an arbitrary distance. To characterize the curvilinear motion of the chains we also define the time $`\tau _{N/2}`$ for the diffusion of the middle monomer over half of the chain’s contour by $`g_6(t\tau _{N/2})=(Nl/2)^2`$. From the analysis it will turn out that two time scales, $`\tau _1`$ and $`\tau _{\text{end}}`$, are particularly important. They have been included in Tables 1 and 2.

#### Finite size effects.

Finally, we note that the results obtained for $`N=1024`$ and $`\varphi =0.5`$ have to be taken with some care. The configuration contains only $`128`$ chains and as $`LR_\text{e}96`$, we cannot rule out dynamical finite size effects — especially for the quenched SSA. Some very large motions have been observed due to chains taking advantage of the image of their own correlation hole in the neighboring periodic boxes. Unfortunately, a systematic study of the influence of the box size $`L`$ is beyond our computational capabilities at present.

## IV Analytical expectations

In the following paragraphes we summarize the relevant theoretical expectations for the analysis of the simulation data.

#### Mean-square displacements: general properties.

If the chains do not become permanently trapped, one expects free curvilinear and spatial diffusion, i.e.,

$`g_6(t)`$ $``$ $`2D_\text{c}t`$ (1)
$`g_1(t)=g_3(t)=g_4(t)`$ $``$ $`6D_\text{s}t,`$ (2)
for times larger than the terminal relaxation time $`\tau _\text{t}`$ and for distances larger than the radius of gyration $`R_\text{g}`$. Furthermore, $`g_2`$ and $`g_5`$, both defined with respect to the center of mass of the chains, must yield a purely static quantity, the radius of gyration,

$`g_2(t)g_5(t)/4`$ $``$ $`R_\text{g}^2,`$ (3)
for very long times in any ergodic dynamical scheme PBHK91\_JCP .

At short times the spatial displacement is directly expressible in terms of the curvilinear displacement $`s`$ along a fractal object characterized by the Flory exponent $`\nu `$. For the MSD of the center of mass this suggests

$$Ng_3(t)N\left(\frac{\stackrel{}{R}_\text{e}}{N}\frac{s(t)}{l}\right)^2\stackrel{\text{MF}}{}\frac{R_\text{e}^2}{N}\frac{g_6(t)}{l^2}(b/l)^2g_6(t)$$
(4)
where the penultimate step follows from a (seemlingly innocent) mean-field assumption. The last step uses the fact that in the melt $`\nu =1/2`$. For the monomer displacements along the backbone a similar reasoning gives

$$g_1(t)=g_2(t)g_4(t)=g_5(t)R(s(t))^2\stackrel{\text{MF}}{}b^2\left(g_6(t)/l^2\right)^\nu $$
(5)
where $`R(s)^2`$ is the size of a piece of chain of curvilinear length $`s`$. We have again made a mean-field approximation and assumed the fractal behavior of asymptotically long chains.

#### Uncorrelated reptation hypothesis.

Until now everything is fairly general and is valid for both the URH and the ARH. For uncorrelated snakes the motion along the backbone must be linear with respect to time. This fixes the short time predictions of the URH:

$`g_6(t)2D_\text{c}t,`$ $`g_3(t)6D_\text{s}t`$ (6)
$`g_1(t)=g_2(t)g_4(t)=g_5(t)`$ $``$ $`t^\nu .`$ (7)
Equation (6) refers, strictly speaking, to the approximately free diffusion along the chain backbone. As our simulation does not deal with phantom chains, there must be a small correction even for very dilute solutions. (An isolated chain is not a phantom chain. It still interacts with itself.) At long times the MSD’s must obey Eqs. (13).

#### Activated reptation hypothesis for a quenched molecular field.

The description of the MSD’s within the ARH is best outlined by considering chains in a molecular field that is quenched at $`t=0`$ (qSSA). Hence, we start with an equilibrium ensemble of reference chains. This implies that the chains initially sit in their correlation hole with relatively low conformational free energies. To creep out of the local, well adapted environment a free energy penalty must be paid on average. It is then reasonable to assume Semenov97 ; SR98  (i) that fluctuations of the initial free energy may be neglected (i.e., we consider a typical initial chain), (ii) that the curvilinear displacement $`s`$ may be taken as the reaction coordinate for the associated Kramers escape problem, and (iii) that the same penalty $`F(s)`$ may be associated with the possibly different conformations characterized by the same $`s`$ — a simple, but less obvious mean-field approximation.

To estimate the dependence of the energy penalty on $`s`$ it was argued in Refs. Semenov97 ; SR98  that $`F(s)=fs^\phi `$ for $`ss_+(N)`$ where $`f`$ is a constant. The upper limit $`s_+(N)`$ depends on the relaxation (escape) mechanism. Since we assumed the molecular field to be quenched, it seems natural to take $`s_+N1`$: At distances larger than the radius of gyration, translational entropy must ultimately win over all other free energy contributions. To close the escape problem we still need to specify the physics for $`s>s_+`$. Because of the gain of the translational entropy $`F(s)`$ should decrease logarithmically at large distances. A crude, but physically sufficient approximation is to relax the penalty completely: $`F(s)0`$ for $`s>s_+`$ Semenov97 .222It is not sufficient to set $`F(s)=F(s_+)`$ for $`s>s_+`$ as this corresponds to a much too high return probability for the one-dimensional Kramers problem into which we have cast the original three-dimensional model.

The exponent $`\phi `$ appears to be a matter of debate. In Refs. Deutsch ; SR98  a linear relationship was suggested ($`\phi =1`$). This implies that the chemical potential for a curvilinear move is constant and thus independent of the tail length $`s`$. On the other hand, it was suggested in Ref. Semenov91 ; Semenov97  that the free energy penalty should result from fluctuations of the molecular field rather than from a fixed chemical potential. The elaboration of this argument yields $`\phi =1/2`$.

Given the free energy barrier one can calculate the mean-square displacements via Eqs. (4,5) if one knows the probability distribution $`p(s,t)`$. In principle, $`p(s,t)`$ can be obtained from the solution of the one-dimensional Smoluchowski equation. Analytically, this calculation has not been done yet,333We determined $`p(s,t)`$ in the simulations, see Fig. 16. but the qualitative behavior of the MSD’s may be inferred from the following arguments. At short times ($`t\tau _{}`$), we have $`ss_+`$, and thus $`F(s)/k_\text{B}T1`$. In this limit, one obviously recovers the free curvilinear diffusion of Eqs. (6) and (7). If time increases, there should be an interval $`\tau _{}t\tau _+`$ such that $`k_\text{B}TF(s)F(s_+)`$. Here, we expect that the chain is temporally localized in its correlation hole, which should lead to subdiffusive behavior due to intermittence of large scale displacements. The subdiffusive regime should extend up to a “hopping time” $`\tau _+`$. $`\tau _+`$ is the time taken by the chain to diffuse over the curvilinear distance $`s_+`$ and thus to overcome the activation barrier. This implies $`\tau _+s_+^2\mathrm{exp}[F(s_+)/k_\text{B}T]`$. Only for $`t\tau _+`$ the chain can diffuse freely.

#### Activated reptation hypothesis for an annealed molecular field.

Up to now, we considered the molecular field to be quenched. If all chains are allowed to move simultaneously, this annealing should reduce the activation barrier. For this case it has been suggested that $`s_+(N)d_{\text{end}}^2N^{2/3}`$ Deutsch ; Semenov97 ; SR98 : The free energy constraints can be relaxed by the reptation of other chains. For short times corresponding to $`ss_+`$ the difference between a quenched or an annealed molecular field should be negligible. In this limit, the correlations between neighboring chains should be essentially of static nature. On the other hand, for $`ss_+`$ the correlations are expected to be due to the cooperative rearrangements of several chains. Semenov and Rubinstein have proposed SR98  an intricate exchange mechanism which allows the relaxation of the free energy penalty between different snakes at $`ss_+`$. It is not obvious to us if one can still cast these processes in the simple form of a one-dimensional Kramers escape problem. In any case this should only affect the probability distributions and MSD’s around $`ss_+`$.

#### Relaxation times and diffusion coefficients.

The relative scaling of the terminal time $`\tau _\text{t}`$ and the diffusion coefficients is imposed by Eq. (4) which fixes the ratio of both diffusion coefficients:

$$\tau _\text{t}\frac{\left(Nl\right)^2}{D_\text{c}}\frac{R_{\text{e}}^{}{}_{}{}^{2}}{D_\text{s}}.$$
(8)
This ensures the matching of the long and short time behavior. Hence, $`\tau _\text{t}(N)`$, $`D_\text{c}(N)`$, $`D_\text{s}(N)`$ are expected to contain the same information. This should hold irrespective of whether the URH or the ARH applies.444In Sec. V.3 it will be checked explicitly that this general relation still holds for the annealed SSA at different densities (Fig. 8). As we shall see, it fails, however, for the quenched SSA at high density (Figs. 11,13).

As we attempt to move all $`N`$ monomers of a snake collectively, the curvilinear diffusion coefficient should be independent of chain length for uncorrelated snakes. Hence, it follows from Eq. (8) and from the URH that

$$D_\text{c}/l^2mN^0\text{ , }D_\text{s}m(R_\text{e}/N)^2\text{ , }\tau _\text{t}\frac{N^2}{m},$$
(9)
$`m`$ being the $`N`$-independent mobility which, due to the URH, should be given by the acceptance rate, $`mAN^0`$. Note that these equations differ by a factor of $`N`$ from the usual ones characteristic of reptating chains within the URH DE86 . This is due to the fact that there, a factor $`N`$ is spent by the internal Rouse dynamics along the primitive path. It is precisely this missing factor which makes the SSA attractive in the first place to explore the dynamics of extremely long chains.

Within the ARH the curvilinear diffusion is fixed by the hopping time $`\tau _+`$ to reach a distance $`s_+`$ within the potential formed by the local correlation hole, i.e.,

$$D_\text{c}(N)/l^2\frac{(s_+/l)^2}{\tau _+}m\mathrm{exp}\left[F\left(s_+(N)\right)\right]=m\mathrm{exp}\left(fs_+(N)^\phi \right)m_+,$$
(10)
which in turn implies $`D_\text{s}(N)`$ and $`\tau _\text{t}(N)`$ according to Eq. (8). Equation (10) is obviously subject to the specific dependences of $`s_+(N)`$ and $`F(s)`$ on their respective parameters. We will try to estimate both $`\tau _+(N)`$ and $`F(s)`$ in Figs. 14 and 16.

## V The dynamics of annealed and quenched snake fields

We begin our presentation of the slithering snake dynamics in Sec. V.1 with the characterization of local phenomenological dynamical scales. Dilute and moderately dense systems are shown in Sec. V.2 to be well described by the predictions of the URH. Diffusion coefficients and characteristic time scales for a broad range of densities are presented in Sec. V.3. We then focus on the dynamical properties of systems of dense solutions and melts. The various MSD’s and their scaling are discussed in the next subsection. By analyzing various displacement histograms possible ergodicity problems are explicitly ruled out for the annealed SSA in Sec. V.5. The histograms of the curvilinear displacements allow us to estimate the free energy penality $`F(s)`$ mentioned above.

### V.1 Local phenomenological dynamical parameters

#### Acceptance rate.

The simplest and statistically most accurate dynamical property sampled in any Monte Carlo simulation is the acceptance rate $`A`$ (= ratio of accepted to proposed moves). It is included in the tables and plotted versus $`N`$ in Fig. 2 and versus $`\varphi `$ in Fig. 3. We note en passant that the acceptance rates of the aSSA and the qSSA are identical. The SSA acceptance rates are also completely independent of chain length. We recall that in contrast a weak, but systematic chain length dependence of $`A`$ is found for local dynamics due to the higher mobility of the chain ends compared to inner monomers (see, e.g., Fig. 10 of Ref. PBHK91  and the bold line in Fig. 2).

As shown in Fig. 3, $`A`$ is density independent below $`\varphi 0.01`$, where it is also much larger for the SSA than for purely local dynamics. This indicates that for these densities the blob size $`\xi `$ becomes sufficiently large to allow dilute SSA dynamics on distances smaller than $`\xi `$. Above $`\varphi 0.1`$ density effects become pronounced and the acceptance rate of the SSA decreases strongly. This decrease may be explained as follows: one monomer is inserted in the SSA, i.e., 8 free lattice sites are needed, whereas 4 empty sites suffice for local moves, in which a monomer is shifted from its original position by one lattice constant. Hence, the observed smaller acceptance rate for the SSA at large density ($`\varphi =0.5`$) is expected.

#### Local mobility.

Obviously, it would be naive to believe that the acceptance rate of the Monte Carlo process necessarily sets the scale for the local mobility of the monomers. There may be many allowed moves not contributing to the motion. (They would also persist in a glass where the mobility is zero.) In order to measure directly the monomeric mobility $`m`$ in the SSA we exploit Eq. (7) and compare the simulation data with

$$g_1(t)=b^2\left(mt\right)^\nu .$$
(11)
This equation provides a reasonable description of the MSD of the central monomer at short times, as Fig. 4 illustrates for a dilute solution and Fig. 5 for a dense melt.

For dilute systems it is appropriate to fit $`g_1(t)`$ using $`\nu =\nu _0`$ and the density independent effective bond length $`b_03`$. The same analysis can also be done in semidilute solutions for times shorter than the time required to cross the blob size so that the excluded volume interaction is not yet screened. This approach was also applied in Ref. PBHK91  to extract $`m`$ for the BFM in the case when only local moves are allowed. Our results for $`m`$ for both the annealed and the quenched SSA are included in Fig. 3.

However, Eq. (11) with $`b=b_0`$ and $`\nu =\nu _0`$ is not the only possible way to do the analysis. For overlapping systems another choice is $`b=b(\varphi )`$ and $`\nu =1/2`$. This measures the effective mobility on distances of several blob sizes. Equation (11) then suggests that $`g_1(t)`$ should scale as $`g_1(t)b(\varphi )^2t^{1/2}`$ at short times, independent of chain length. Figure 5 confirms this expectation. If the mobility is now extracted from Fig. 5, the differences between the new values and those obtained before (using $`b=b_0`$ and $`\nu =\nu _0`$) are negligible.

Here, the important point is that the local mobility is the same for all chain lengths at a given density. (Note that the mobilities obtained by Paul et al. PBHK91  for the local dynamics are also independent of chain length.) Furthermore, Figs. 4 and 5 also show that $`m`$ is independent of whether the system is annealed (aSSA) or quenched (qSSA). For both dynamical schemes, the mobility rapidly decreases with density if $`\varphi >0.1`$. This behavior is qualitatively the same as that of the acceptance rate. Similar results are obtained if Eq. (11) is fitted to the MSD of the chain ends $`g_4(t)`$ (not shown). For the SSA, the chain ends are more mobile than the central monomer by a factor of order one.

In summary, the acceptance rate or the local mobility cannot explain the unexpected chain length dependence of the motion of strongly entangled chains to be discussed in Secs. V.3V.5.

#### Effective mobility for the aSSA.

A glance at Fig. 5 shows that in fact two fits using Eq. (11) are possible for the annealed SSA once $`\varphi \varphi ^{}`$. The rescaling $`y=g_1(t)/b(\varphi )^2t^{\nu =1/2}`$ yields a data collapse at very short times, $`t\tau _{}`$. From the resulting plateau the local mobility $`m`$, shown in Fig. 3, was extracted. For larger times, however, $`y`$ decreases within the window $`\tau _{}t\tau _+(N)`$. The existence of such a time window is obviously at variance with the uncorrelated reptation hypothesis. For still larger times, $`\tau _+t\tau _\text{t}`$, the rescaled MSD becomes approximatively horizontal. This allows us to define an effective mobility $`m_+(N)`$ (included in the tables). It reflects the complicated correlations in the second regime $`\tau _{}t\tau _+`$, which are responsible for its pronounced dependence on $`N`$. The importance of $`m_+`$ for the description of diffusion coefficients and relaxation times will be clarified in Sec. V.3. Finally for $`t>\tau _\text{t}`$, Fig. 5 shows that $`y\sqrt{t}`$, implying that the monomers diffuse freely.

### V.2 Dilute and weakly overlapping systems

Before analyzing the intricate physics at high chain overlap we briefly consider the simpler situation of dilute or weakly overlapping solutions. These systems serve as a reference point for the following discussion of concentrated solutions and melts. That the predictions of the URH apply in the limit of high dilution can be inferred from Figs. 6, 7 and 9. The rescaling used for the vertical axes of the three figures is motivated by Eq. (9) and the fact that the local mobility $`mA(\varphi )N^0`$. The figures confirm the validity of the URH by showing that $`D_\text{c}N^0`$, $`ND_\text{s}N^{2\nu _01}`$ and $`\tau _\text{t}\tau _1N^2`$.

Figure 4 presents the spatial MSD’s for a typical dilute or weakly overlapping configuration. We compare the measured MSD’s with the expected long time behavior (Eqs. (2,3)) and the short time predictions of the URH (Eqs. (6,7) using $`\nu =\nu _0`$). Note that the MSD of the center of mass $`g_3(t)`$ is linear over several orders of magnitude, implying that corrections due to excluded volume correlations of the snake with itself are very small. The same applies to the curvilinear motion $`g_6(t)`$ (not shown). Incidentally, the excellent agreement confirms that our program is running properly and that the URH is appropriate as long as chain interactions are weak. The presented data are for the annealed SSA, but absolutely identical results have been obtained for qSSA, as expected.

The above statements remain valid for all configurations and both dynamical schemes provided that $`\varphi \varphi ^{}`$.

### V.3 Diffusion constants and time scales

#### Diffusion coefficients for the aSSA.

The curvilinear and spatial diffusion coefficients, $`D_\text{c}`$ and $`ND_\text{s}`$, are displayed in Figs. 6 and 7 not only for dilute solutions, but also for densities ranging up to the melt density $`\varphi =0.5`$. We plotted $`D_\text{c}/A`$ and $`D_\text{s}N/Ab^2`$ versus chain length, since we expected to find more or less horizontal curves (see Eq. (9)).

Surprisingly, Eq. (9) breaks down with increasing chain overlap. For a given density the data points become progressively curved in log–log coordinates. In fact, $`D_\text{c}`$ and $`ND_\text{s}`$ drop over more than two decades at $`\varphi =0.5`$, as $`N`$ increases. It is important to stress that the chain length dependence of $`D_\text{c}`$ and $`ND_\text{s}`$ increases gradually with chain overlap and that this effect is by no means a pathology of a lattice model at high density.

The data for $`\varphi =0.5`$ may be compared to a phenomenological power law with exponent $`\beta =4/3`$. This power law is merely a guide to the eye and a rough attempt to characterize the slope in log–log coordinates for the available chain lengths. However, the stretched exponential $`D_\text{c}ND_\text{s}\mathrm{exp}(0.8N^{1/3})`$ motivated by Eq. (10) describes the data much better. Our scaling is, in fact, in excellent agreement with the ARH and the proposed energy barrier exponent $`\phi =1/2`$ Semenov97 . The alternatively suggested value $`\phi =1`$ Deutsch ; SR98  is not compatible with the simulation data (at least for the chain lengths simulated up to now).

If, instead of the acceptance rate, we use $`m_+(N)`$ to rescale the diffusion coefficient, i.e., $`ND_\text{s}/m_+b(\varphi )^2`$, we find a very good data collapse (small symbols at the top of Fig. 7). A similar collapse is also possible for $`D_\text{c}/m_+`$ (not shown). This shows that it is possible to absorb the complicated physics at short times into one parameter, $`m_+(N)`$. Thus, the dynamics for $`\tau _+t\tau _\text{t}`$ is perfectly described in terms of uncorrelated reptation dynamics along the fractal chains. It remains to describe and possibly explain the dynamics in the second regime, $`\tau _{}t\tau _+`$, of Fig. 5.

#### Diffusion coefficients for the qSSA.

With increasing chain overlap the qSSA diffusion coefficients decrease even more rapidly with density and chain length than their aSSA counterparts. In fact, while we have obtained most of the curvilinear diffusion coefficients from our qSSA simulations, we have been unable to measure $`D_\text{s}`$ for $`\varphi >0.125`$. This surprising statement is obviously at variance with Eq. (8). It indicates that curvilinear and spatial motions decouple for these high densities due to the localization of the snakes within their correlation hole in the quenched case. Furthermore, at high density $`D_\text{c}`$ is found to be a power law without any detectable curvature (Fig. 6). However, it is interesting to check whether Eq. (10) with $`s_+(N)=N`$ for the qSSA and $`\phi =1/2`$ (for both schemes) yields a reasonable fit. This is indeed the case: $`D_\text{c}/A20\mathrm{exp}(0.15N^{1/2})`$ is hardly distinguishable from the indicated slope $`\beta 1.5`$. Longer chains are obviously necessary to confirm the exponential behavior expected by the ARH.

#### Relaxation times for the aSSA.

In Figs. 8 and 9 we display some of the characteristic times defined at the end of Sec. III. In the first figure we test Eq. (8) by plotting various ratios $`\tau _\text{i}/\tau _1`$. Apparently, the ratios are independent of chain length. We find $`\tau _\text{s}/4\tau _\text{c}/5\tau _1\tau _2/0.85\tau _3/4.7\tau _4/0.3\tau _{N/2}/6`$. Only the data for $`\varphi =0.5`$ have been included, but the same results prevail for all densities. This demonstrates that all indicated times are proportional to one characteristic time scale, the “terminal time” $`\tau _\text{t}`$, containing the same information as the diffusion coefficients. That is, Eq. (8) still holds, although Eq. (9) fails.

Accordingly, Fig. 9 shows, e.g., that $`\tau _1`$ is proportional to $`N^2`$ in the dilute limit, but increases like $`\tau _1N^2/D_\text{c}N^2\mathrm{exp}(0.8N^{1/3})`$ in the melt. The figure also includes $`\tau _d`$ for $`d=3`$ and $`\tau _{\text{end}}`$. The first time indicates the distance $`d`$ up to which $`g_1(t)`$ remains independent of $`N`$, and characterizes the extent of the first dynamical regime displayed in Figs. 5 and 12. In fact, it strongly decreases with monomer density, and we expect it to scale like the time needed to diffuse over the excluded volume blob $`\xi `$. From the observed behavior of $`D_\text{c}`$ one anticipates for $`\tau _{\text{end}}`$ the scaling $`\tau _\text{t}\tau _{\text{end}}s_{\text{end}}^2/D_\text{c}N^{2/3\nu }\mathrm{exp}(0.8N^{1/3})`$, where $`s_{\text{end}}d_{\text{end}}^{1/\nu }`$. This is again consistent with the data presented in the figure.

#### Relaxations times for the qSSA.

Data points for $`\tau _1`$ and $`\tau _{N/2}`$ at $`\varphi =0.5`$ have been included for the qSSA in Fig. 9. Unfortunately, the range of $`N`$, which could be simulated, is insufficient to check if indeed $`\tau _1N^2\mathrm{exp}(\text{const}N^{1/2})`$. The time scale $`\tau _{N/2}`$, measuring the curvilinear motion over $`N/2`$ monomers, shows an effective power law which is consistent with the corresponding diffusion coefficient from Fig. 6: $`\tau _{N/2}N^2/D_\text{c}N^{2+1.5}`$. Note that $`\tau _1`$ and $`\tau _{N/2}`$ scale differently in contrast to what was observed for the annealed SSA (Fig. 8).

### V.4 Mean-square displacements of dense snakes

#### Long time limit.

Focusing on systems in the melt limit ($`\varphi =0.5`$) we now discuss several MSD’s and their scaling in detail (Figs. 1014). As a glance at the presented MSD’s shows, the long time predictions (Eq. (13)) are in very good agreement with the data for the aSSA, indicating that the chains do not get trapped and that the simulation runs have been long enough. This can be seen in Fig. 10 for the spatial MSD’s for $`N=512`$, in Fig. 11 for both spatial and curvilinear MSD’s for $`N=128`$ and in Fig. 12 for the curvilinear MSD $`g_6(t)`$.

The situation is more complex for the data of the qSSA (included in Figs. 1112). Figure 11 shows that the spatial MSD’s increase (approximately) logarithmically for the longest times we have been able to simulate. This must be a transient: For very long times these MSD’s must either become constant, i.e., the chains become localized on a length scale (presumably) given by the radius of gyration, or ultimately, they diffuse freely in space. So far, we have been able to show conclusively only the correctness of the second proposition for small chains ($`N32`$).

At first sight surprisingly, it is much less difficult to reach the free curvilinear diffusion limit, as may be seen from Fig. 11, where $`g_6(t)`$ and $`Ng_3(t)`$ are directly compared for $`N=128`$, and from Fig. 12, where the $`g_6(t)`$ for different masses are presented. This indicates that for the qSSA the curvilinear and spatial MSD’s decouple for $`t\tau _{\text{coup}}`$ and $`ss_{\text{coup}}`$ (see Fig. 11); although unable to move much further in space, the snakes manage to move relatively rapidly within the cage of their correlation hole. This is why we prefer to call this surprising behavior of dense qSSA snakes (permanent or transient) localization in contrast to a trivial local blocking of the snake moves (on a scale given by the monomer density), which must obviously occur at very high densities.

#### Short time limit for the aSSA.

While the MSD’s are well described by the power laws of the URH in the dilute limit (Fig. 4), this becomes strikingly different at larger densities in the time window $`\tau _{}t\tau _+`$. This is revealed by the direct comparison of dilute and dense systems with $`N=512`$, as given in Figs. 4 and 10. In this paper we attempt to characterize and to explain the features visible in the latter plot.

We have tentatively described the subdiffusive behavior of the aSSA data by means of an effective power law:

$$Ng_3(t)g_1(t)^2g_6(t)t^\alpha \text{ with }\alpha 1/2.$$
(12)
This power law provides a remarkable fit for $`\tau _{}t\tau _+`$ (compare Figs. 5,10,12), even though some curvature is evident in log–log coordinates. The alternative fit

$$Ng_3(t)g_1(t)^2g_6(t)\mathrm{log}(t)^{4\alpha }$$
(13)
is less satisfactory for our aSSA data (fitting a smaller time interval), as may be seen from Fig. 11. (The exponents in the two previous equations have naturally been chosen consistently with the general relationship expressed by Eqs. (4,5) using $`\nu =1/2`$.)

#### Short time limit for the qSSA.

Equation (13), however, reasonably fits the qSSA data. This may be seen from Figs. 5 and 11. While for the spatial MSD’s this fit works for all times $`t\tau _{}`$, $`g_6(t)`$ becomes linear for $`t\tau _{+\text{q}}`$, $`t\tau _{+\text{a}}`$ with $`\tau _{+\text{q}}\tau _{+\text{a}}`$ (compare aSSA and qSSA in Fig. 12 for $`N=64,128`$ and $`1024`$), where we altered the index to distinguish the different algorithms. Note that annealed SSA and quenched SSA data are identical for small times $`t\tau _{\text{fluc}}(N)`$ (Figs. 11,12). This is due to the low probability that a snake in the aSSA may interact with the density fluctuations created by neighboring snakes at short times. Hence, the molecular fields probed behave in the same way. As $`\tau _{\text{fluc}}(N)`$ increases strongly with $`N`$ (see below), this suggests that the logarithmic behavior might be seen for aSSA snakes as well, if we could simulate longer chains.

#### Scaling of $`g_1(t)`$, $`g_3(t)`$ versus $`g_6(t)`$.

In Fig. 13 we verify now directly whether the time dependence of all spatial MSD’s is indeed fully encapsulated in terms of the curvilinear MSD $`g_6(t)`$, as suggested by both the URH and the ARH. In the main panel of the figure we check this for the aSSA data at $`\varphi =0.5`$ by plotting $`g_1(t)`$ and $`g_3(t)`$ versus $`g_6(t)`$. The MSD of the center of mass is linear with regard to $`g_6`$, as expected. This is in accord with the previous finding $`D_\text{c}ND_\text{s}`$ for Gaussian chains. We find a similar data collapse for $`g_1(g_6)`$ for $`t\tau _\text{t}`$. For sufficiently long chains, we clearly recover the anticipated power law of Eq. (7) for the diffusion on an object of inverse fractal dimension $`\nu =1/2`$ (dashed line in Fig. 13). (The departure from the slope $`1/2`$ for small $`g_6(t)`$ is due to the non-gaussian statistics for small distances along the chain.) We conclude that the curvilinear and spatial diffusion are translated into each other, as expected, and that the departure from the URH can be traced back to the anomalous curvilinear diffusion depicted in Fig. 12.

In the inset of Fig. 13, we present $`Ng_3(g_6)`$ for the qSSA for two chain lengths and two densities. While for the lower density $`g_3(t)`$ remains a (roughly) linear function of $`g_6(t)`$, this is only true for short times at $`\varphi =0.5`$, where the spatial displacements increase much more strongly. Note that both densities appear to be independent of chain length. This indicates a very weak mass dependence for $`s_{\text{coup}}`$.

#### Scaling with chain length.

We now discuss in more detail the chain length dependence of the subdiffusive behavior at fixed volume fraction $`\varphi =0.5`$. Since we have shown (Fig. 13) that the spatial MSD’s of the annealed SSA scale in terms of $`g_6(t)`$, we can concentrate on the $`N`$ dependence of the latter quantity (Fig. 12). The vertical axis is normalized to yield the curvilinear diffusion coefficient $`D_\text{c}`$ as the plateau value for long times. At short times $`t\tau _{}50`$, all curves merge. Hence, $`\tau _{}`$ is chain length independent, but depends weakly on the density, as may be shown by a similar plot containing data for different $`\varphi `$. For times beyond $`\tau _{}`$ the dynamics slows down more and more until it reaches the plateau of free diffusion at $`\tau _{+\text{a}}(N)`$ for aSSA and $`\tau _{+\text{q}}(N)`$ for qSSA, respectively. Apparently, the correlations become more and more pronounced with increasing $`N`$.

We have directly computed the four times $`\tau _{+\text{a}}`$, $`\tau _{+\text{q}}`$, $`\tau _{\text{fluc}}`$ and $`\tau _{\text{coup}}`$ and the related curvilinear displacements. Our measurements are not very precise (particularly, for the qSSA data better statistics is needed), and we only present some general trends here. A comparison of all chain lengths shows that

$`\tau _{\text{fluc}}(N)`$ $``$ $`\tau _{+\text{a}}(N)\tau _{\text{end}}(N)`$ (14)
$`\tau _{\text{coup}}(N)`$ $``$ $`\tau _{+\text{q}}(N)N^2.`$ (15)
In other words, both pairs of time scales correspond to different definitions of the same physical processes. In the first case, these are the density fluctuations in the annealed scheme, causing aSSA and qSSA to differ and the curvilinear motion to become eventually diffusive. The second equation corresponds to the free curvilinear diffusion of the snake trapped in its quenched correlation hole. The power law indicated in Eq. (15) is an estimate motivated by the data given in Fig. 9 (diamonds).

The last claim of Eq. (14), which is central to the understanding of the observed correlations and is in line with the ARH Deutsch ; SR98 , is further tested in Fig. 14. The attempt $`\tau _{+\text{a}}\tau _1`$ fails to scale the crossover regime between the subdiffusive time window and the free curvilinear diffusion at late times. (Note that there is always the time scale $`\tau _{}N^0`$ for the free diffusion at very short times and no scaling covering all regimes is possible.) Instead, the alternative $`\tau _{+\text{a}}\tau _{\text{end}}`$ is successful. The observed scaling confirms Eq. (14) and demonstrates that the crossover from correlated to mean-field behavior is determined by $`d_{\text{end}}`$.

Finally, we note for the second, less important, relaxation process for the qSSA that $`s_{+\text{q}}`$ increases very weakly with mass, $`s_{+\text{q}}^2N^{0.25}`$. (However, this does not contradict the scaling documented in the inset of Fig. 13. Note also that this is again a rough estimate, but we definitely have to rule out the two “natural” guesses $`s_{+\text{q}}N`$ or $`s_{+\text{q}}N^0`$.) Interestingly, $`s_{+\text{q}}^2N^{0.25}`$ and $`\tau _{+\text{q}}(N)N^2`$ are compatible with the (very precisely measured) qSSA diffusion coefficient from Fig. 6, $`D_\text{c}s_{+\text{q}}^2/\tau _{+\text{q}}1/N^{1.5}`$, and, hence, with the time scale $`\tau _{N/2}N^{3.5}`$ shown in Fig. 9. While this is satisfactory, we presently cannot explain the scaling behavior of $`\tau _{+\text{q}}`$ and $`s_{+\text{q}}`$.

### V.5 Histograms and Non-Gaussianity

So far, we have discussed results obtained from the sampled mean-square displacements and related diffusion constants and time scales. We now turn to the additional information contained in the histograms and higher moments of the curvilinear displacements. Merely results from the annealed field dynamics are discussed here, since statistics is too poor for the qSSA.

An important, not only technical, question concerns the possibility that some chains or monomers become locally completely blocked. We have paid much attention to this issue and have carefully checked that, in fact, all chains and monomers ultimately move arbitrarily far in space for the densities studied ($`\varphi 0.5`$). Clearly, for some higher density this must eventually break down.

One may readily construct by hand chain configurations of the BFM, which lead to blocked chains for any dynamical scheme under consideration, even if local moves are allowed. These extremely rare points of the phase space also exist for our slithering snake scheme, although they are apparently excluded from our initial configurations (which is the physically important point to make here). These states cannot be accessed (by construction) and in this strict mathematical sense our dynamics is not ergodic. However, this is completely irrelevant for the thermodynamic and dynamic behavior.555 In this context it is also of relevance that one obtains the same dynamics if a small amount of local hopping moves, $`\omega 1`$, is included in the algorithm which allows additional transitions between points of the phase space. See Fig. 18.

#### Histogram of first passage time.

Two checks are presented in the histograms given in Figs. 15 and 16. The first figure presents the histograms of first passage times of the central monomer over a curvilinear distance $`s=lN/2`$ for various chain lengths and densities. At long enough times, typically $`10\tau _{\text{fp}}`$, where $`\tau _{\text{fp}}=t_{\text{fp}}`$ is the mean first passage time, all chains have crept out of their initial tube. If properly rescaled, the data collapse onto a master curve. The distribution of the first passage times can be calculated if free curvilinear diffusion is assumed. It is a linear superposition of modes $`p`$ with an exponential time dependence $`\mathrm{exp}(p^2t_{\text{fp}}/\tau _\text{c})`$. For long times, this sum is dominated by the lowest mode and, accordingly, we find an exponential decrease $`\mathrm{exp}(t_{\text{fp}}/\tau _\text{c})`$ for all systems. As expected, the mean first passage time scales like the terminal time, $`\tau _{\text{fp}}6\tau _15\tau _\text{c}/6`$, and is of the same order as $`\tau _{N/2}`$ which was defined as the time where $`g_6(t)`$ corresponds to the same curvilinear distance.

#### Histogram of curvilinear displacements $`s`$.

Figure 16 presents the histogram of $`s(\mathrm{\Delta }t)`$ for $`N=1024`$ and $`\varphi =0.5`$ for various time intervals, as indicated. We focus on the interesting time window $`\tau _{}t\tau _+\tau _\text{t}`$. First of all we stress that the histograms have (very qualitatively) the same shape and fall off monotonously without any singularity at the origin or elsewhere, i.e., no snake gets blocked at small $`s`$. If no correlations were present, the rescaled histograms should all collapse onto a single Gaussian. More importantly, however, Fig. 16 shows that the different dynamical regimes correspond to (slightly) different distributions. While at times $`t\tau _{}`$ and $`t\tau _+(N=1024)0.210^8`$ the distributions are in fact Gaussian666For clarity, both regimes are not included in the figure. Note that Gaussian behavior is only found in the first regime at much smaller densities than the one discussed in Fig. 16. A sufficiently large blob is required to obtain a Gaussian distribution., they are essentially exponential in the subdiffusive regime where density fluctuations are important. The distribution is thus much broader than originally expected, but all moments of the distribution exist. Although this is not an equilibrium distribution measuring directly the free energy as a function of $`s`$, the Smoluchowski equation suggests that $`\mathrm{log}(p(s))`$ should be proportional $`F(s)`$. We have therefore compared the two proposed exponents for the activation penality $`\phi =1`$ (bold line) and $`\phi =1/2`$ (dots) with the data. Unfortunately at large $`s`$, the statistics becomes too poor to distinguish between both exponents. More work is required to relate the form and scaling of the observed histograms to the respective MSD.

#### Non-Gaussian parameter.

The non-gaussian parameter $`\alpha _2\frac{1}{3}\frac{s^4}{s^2^2}1`$ gives a compact means to describe the Gaussianity of these distributions as a function of time. This is presented in Fig. 17. We recall that $`\alpha _2=0`$ for a Gaussian and $`\alpha _2=1`$ for an exponential distribution. At short times $`\alpha _2(t)`$ decreases rapidly and is chain length independent. The non-gaussianity at these times is due to local monomer interactions. In agreement with the histograms shown in Fig. 16, $`\alpha _2(t)`$ increases strongly for larger times. This shows again that correlations become more pronounced. The effect increases strongly with chain length in agreement with the subdiffusive behavior in Fig. 12. Note that $`\alpha _2`$ becomes even much larger than $`1`$ for our largest chain $`N=1024`$, reflecting the existence of a broad tail in the distribution shown in Fig. 16. We caution, however, that finite size effects might contribute to this effect, as the end-to-end distance of this configuration is comparable to the box-size. It is possible that some snakes manage to find their way rapidly through the periodic box by taking advantage of their own groove spanning the box. For times larger than $`\tau _+(N)`$ the non-gaussianity vanishes again for all $`N`$, as it should.

## VI Summary and outlook: snake versus local moves

#### Summary.

In this paper we reported first numerical results on the dynamics of monodisperse polymer chains moving according to the slithering snake algorithm SSA (no local relaxation pathways). The simulations were done with the three-dimensional bond-fluctuation model (BFM). We investigated two different dynamical schemes by focusing mainly on the polymer melts (high monomer densities): Either we allow all snakes to move (aSSA), which permits fluctuations of the molecular field at the ends of neighboring chains, or all snakes but one are quenched (qSSA). This work concentrates on the annealed algorithm, whereas the second mainly serves as a benchmark to investigate the role of density fluctuations due to neighboring chains. We studied the scaling of various spatial and curvilinear mean-square displacements (MSD), diffusion coefficients and relaxation times with chain length. Furthermore, we also presented the distributions of the curvilinear motion.

Our results are broadly in agreement with the “activated reptation” hypothesis (ARH) Deutsch ; SR98 . For the aSSA, we find that at high chain overlap the terminal relaxation time $`\tau _\text{t}`$ increases much more strongly with chain length $`N`$ than possible for uncorrelated snakes which are obtained for dilute and weakly overlapping chains. In fact, a stretched exponential $`\tau _\text{t}N^2\mathrm{exp}(0.8N^{1/3})`$ fits our data well. Similar and directly related (Eq. (8)) behavior can be found for the spatial and curvilinear diffusion coefficients, $`ND_\text{s}D_\text{c}`$, and other time scales, such as $`\tau _{\text{end}}`$, which measures the time needed for the diffusion over the mean distance between chain ends $`d_{\text{end}}`$. The stretching exponent $`1/3`$ is in agreement with the free energy penalty $`F(s)s^{1/2}`$ ($`s`$ being the curvilinear displacement) estimated in Ref. Semenov97 , but is incompatible with the more recent suggestion $`F(s)s`$ of Ref. SR98  (at least for the chain lengths accessible in the present simulation).

A more detailed study of different moments of the displacements shows that the snake dynamics is only uncorrelated (exhibiting a free curvilinear diffusion with Gaussian distribution) for times smaller than $`\tau _{}(\varphi )N^0`$, which is determined by the blob size $`\xi `$, and for times larger than $`\tau _+\tau _{\text{end}}`$ (shown by the scaling collaps of Fig. 14), which is due to the density fluctuations present for distances larger than $`d_{\text{end}}`$. They are responsible for an additional time scale (other than the terminal time) in the annealed algorithm, allowing the chains to relax the free energy penality caused by their motion in a locally rigid and frozen network of obstacles retaining them in their correlation hole. This is supported by the finding that the MSD’s of both the aSSA and the qSSA, which are strictly identical at small times, become different around $`\tau _+`$.

In the intermediate time interval the dynamics becomes increasingly subdiffusive and non-gaussian (Figs. 16,17), as the chain length is increased for both aSSA and qSSA dynamics. The precise description of the MSD’s in this regime is not evident for the masses we have been able to simulate and our study is not fully conclusive. While power laws yield nominally better fits for aSSA, logarithmically slow relaxation processes might also be appropriate. The observed behavior for the quenched dynamics (Figs. 11,12), which could provide an “envelope” that the annealed dynamics might follow more closely with increasing chain length, is indeed better described in this way. Logarithmic behavior is, however, not easily reconciled with (mean-field) ARH approach which predicts an algebraic relationship. More simulations are warranted to clarify this issue.

The curvilinear displacement histograms decrease monotonously without any singularity at small $`s`$. In the subdiffusive intermediate time they decrease essentially exponentially (Fig. 16). This suggests that the free energy penalty $`F(s)`$ is indeed of power law form. A direct measurement of the acceptance rate to perform a step $`ss+1`$ is currently underway using the quenched scheme which should allow a direct computation of the free energy cost for an escape attempt of order $`s`$.

In summary, we believe that our study demonstrates that the activated reptation approach is relevant for the description of concentrated solutions and melts of purely slithering snakes without local motion ($`\omega =0`$).

#### The mapping onto a locally realistic model.

It is important to recall that the slithering snake dynamics discussed in this paper is highly artificial, as no local rearrangements are allowed and density fluctuations can only occur at chain ends. The idea is that the snake corresponds to the primitive path of the reptation model and each monomer to a coarse-grained tube segment of an originally local model. The mapping of both levels of description is non-trivial. It is not evident whether it can be performed at all without an additional tuning parameter, such as the frequency of local hopping moves.

The mapping might involve the density. It is possible that some of the high molecular densities, which we studied here, are in fact too large to correspond to realistic volume densities of an underlying microscopic model. A characteristic density might therefore exist below which even asymptotically long snakes remain uncorrelated. At present, our computational results do not seem to favor this option, as suggested, e.g., by Fig. 7. However, it remains a possibility which is currently pursued using lower densities and much larger chains such that $`1\varphi \varphi ^{}(N)`$. In any case, a more detailed analysis of the density dependence is warranted to establish the scaling of $`\tau _{}(\varphi )`$ and $`\tau _+(N,\varphi )`$. In other words, it has to be confirmed unambiguously that the first time scale is related to the blob size $`\xi (\varphi )`$ and the second one to the mean distance between chain ends $`d_{\text{end}}(N,\varphi )`$.

#### On the importance of local motion.

Another issue is related to the fact that our slithering snakes do not have any mechanism to relax constraints via local pathways. It might well be that some local dynamical degrees of freedom are necessary to tune a mesoscopic computation scheme such that it becomes consistent with real reptational dynamics. A finite fraction of local moves should weaken the strength of the obstacles imposed by neighboring chains. This is indeed borne out by preliminary studies presented in Fig. 18, where the rescaled spatial diffusion coefficient is plotted for various ratios of local to snake moves $`\omega `$. For short chains ($`N<64`$) $`ND_\text{s}`$ decreases monotonously with increasing $`\omega `$, since local moves are less efficient in exploring the phase space and the confinement is neglible. As the mass increases, we find a striking non-monotonous behavior. The dynamics first becomes more rapid, as local dense fluctuations become more relevant. This effect apparently saturates at a frequency of order $`10`$, which causes the diffusion coefficient to decrease again strongly with increasing $`\omega `$ (at fixed $`N`$). A finite $`\omega `$ generates a lateral tube of size $`d_\text{l}(\omega ,\varphi )>a`$, which should become identical to the tube diameter $`d_\text{e}(\varphi )`$ of the reptation model in the limit of very large $`\omega `$. Recalling that $`d_\text{e}(\varphi )\xi (\varphi )`$ we speculate that a successful mapping of a SSA scheme on a local scheme may require the tuning of the ratio $`d_\text{l}(\omega ,\varphi )/\xi (\varphi )=\text{const}1`$, which implies a different $`\omega `$ for each density. We are currently pursuing further simulations to clarify these questions.

#### A final word.

It might well turn out that even in this larger modeling space no $`\omega (\varphi )`$ can be found, where $`D_\text{c}ND_\text{s}`$ remains chain length independent for large $`N`$. We would regard this as a strong evidence for the correctness of the ARH for a model with fully realistic local dynamics. We believe, however, that astronomically long chains are required to obtain measurable effects.

###### Acknowledgements.

During the course of this work, we had valuable discussions and obtained insights from many colleagues, especially from M. Müller, A. Johner, M. Fuchs, A. N. Semenov and S. Obukhov. LM and JPW acknowledge computational support by IDRIS/France. EL acknowledges support through a fellowship from the Max Planck Institute for Polymer Research (Mainz, Germany).