Reptation in the Rubinstein-Duke model: the influence of end-reptons dynamics

12 Pages·0.35 MB·English

by  Enrico Carlon

Reptation in the Rubinstein-Duke model: the influence of end-reptons dynamics
Enrico Carlon1, Andrzej Drzewin´ski2 and J.M.J. van Leeuwen3
1Theoretische Physik, Universit¨at des Saarlandes, D-66123 Saarbru¨cken, Germany
2Institute ofLow Temperature andStructure Research, PolishAcademy of Sciences, P.O.Box 1410, 50-950 Wrocl aw 2, Poland
2 3Instituut-Lorentz, University of Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands
0 (February 1, 2008)
0
We investigate the Rubinstein-Duke model for polymer reptation by means of density-matrix
2
renormalizationgrouptechniquesbothinabsenceandpresenceofadrivingfield. Intheformercase
n
the renewal time τ and the diffusion coefficient D are calculated for chains up to N =150 reptons
a and their scaling behavior in N is analyzed. Both quantities scale as powers of N: τ ∼ Nz and
J
x
D∼1/N withtheasymptoticexponentsz =3andx=2,inagreementwiththereptationtheory.
1
For an intermediate range of lengths, however, the data are well-fitted by some effective exponents
3
whose values are quite sensitive to the dynamics of the end reptons. We find 2.7 < z < 3.3 and

1.8 < x < 2.1 for the range of parameters considered and we suggest how to influence the end
]
h reptons dynamics in order to bring out such a behavior. At finite and not too small driving field,
c weobservetheonset oftheso-called bandinversion phenomenonaccording towhichlongpolymers
e
migrate faster than shorter ones as opposed to thesmall field dynamics. For chains in the range of
m
20reptonswepresentdetailed shapesof thereptatingchainasfunctionofthedrivingfieldandthe
- end repton dynamics.
t
a
t PACS numbers: 83.10.Nn, 05.10.-a, 83.20.Fk
s
.
t
a
m I. INTRODUCTION since it is a technique which allows to separate polymers
according to their length.
-
d The previous examples demonstrate how reptation is
The idea of reptation was introduced about thirty
n an important mechanism for polymer dynamics in dif-
years ago by de Gennes [1] in order to explain some dy-
o
ferent physical situations. A prominent role in under-
namical properties of polymer melts of high molecular
c
standing the subtle details of the dynamics was played
[ weight. The dynamics of such systems is strongly influ-
by models defined on the lattice, which besides describ-
enced by entanglement effects between the long polymer
1 ing correctly the process at the microscopic level, offer
chains. The basic idea of reptation is that each polymer
v
important computational advantages. The aim of this
3 is constrained to move within a topological tube due to
paper is to investigate in detail one of such lattice mod-
8 the presenceofthe confining surroundingpolymers[2,3].
els, which was originally introduced by Rubinstein [11]
5 Within this tube the polymer performs a snake-like mo-
1 tion and advances in the melt through the diffusion of and later extended by Duke [12] in order to include the
0 stored length along its own contour. One focuses thus effectofanexternaldrivingfield. TheRubinstein-Duke
2 (RD) model has been studied in the past using several
on the motion of a single test chain, while the rest of
0 techniques; there are a limited numbers of exact results
the environmentis consideredfrozen,i.e. as formedby a
/
t networkoffixedobstacles. Thereptationtheorypredicts available [13–15] while a richliterature onsimulation re-
a
sults on the model exists. The latter have been mostly
m thattheviscosityµandlongestrelaxationtimeτ (known
as renewal time) scale as µ τ N3, where N is the obtainedbyMonteCarlo(MC)simulations[11,12,16,17].
- ∼ ∼ MC simulations are tedious for long chains since the re-
d length of the chains, while the diffusion constant scales
n as D 1/N2. These results are not far from the exper- newal time scales as N3, thus it is difficult to obtain a
o iment∼al findings. Measurements of viscosity of concen- small statistical error for large N.
c trated polymer solutions and melts of different chemical In this paper we study the RD model by means of
:
v composition and nature are all consistent with a scal- Density Matrix Renormalization Group (DMRG) [18],
i ing µ N3.4 [4], while for the diffusion constant both a technique which has been quite successful in recent
X
D 1/∼N2 [5] and D 1/N2.4 [6–9] have been reported. years, in particular in application to condensed matter
r ∼ ∼
a This discrepancy triggered a substantial effort in order problems as quantum spin chains and low-dimensional
to reconcile theory and experiments. strongly correlated systems [19]. Using the formal sim-
Another physical situation where reptation occurs is ilarity between the Master equation for reptation and
in gel electrophoresis, where charged polymers diffuse the Schr¨odinger equation, DMRG also allows to calcu-
throughtheporesofagelundertheinfluenceofadriving lateaccuratelystationarystatepropertiesoflongreptat-
electricfield[10]. Thegelparticlesformafrozennetwork ing chains.
ofobstaclesinwhichthepolymermovesthroughthedif- Some of the results reported here, in particular con-
fusion of stored length. Electrophoresis has important cerning the scaling of the renewaltime τ, havebeen pre-
practical applications, for instance in DNA sequencing, sented before [20]. Here we will give full account of the
1
details of the calculations and presenta series of new re- Whentwoormorereptonsaccumulateatthesamelat-
sults concerning reptation in the presence of an electric ticesitetheyformapartofstoredlength,whichcanthen
field. One of the main conclusions of our investigationis diffuse along the chain. In terms of relative coordinates
that the exponents describing the scaling of τ and D in a segment of stored length corresponds to y = 0, there-
i
the intermediate length region appear to be rather sen- fore allowed moves are interchanges of 0’s and 1’s, i.e.
sitive on the structure of the end repton. We find that, 0, 1 1,0. On the contrary, the end reptons of the
± ↔ ±
byinfluencingtheendreptondynamics,onehasaregime chaincanstretch(0 1)orretracttothesiteoccupied
→±
wheretheeffectiveexponentsareconsiderablylowerthan by the neighboring repton ( 1 0). The dynamics of
± →
the standardly found experimental values. Experiments thechainisfullyspecifiedoncetheratesforthemovesare
aresuggested,involvingpolymerarchitecturesnotinves- given. A field ε along the projection axes is introduced
tigated so far, and which ought to confirm our predic- so that moves forward and backward in the field occur
tions. with rates B = exp(ε/2) and B−1 = exp( ε/2). In the
−
This paper is organizedas follows: in Sec. II we intro- following we will be interested in both cases ε = 0 and
duce the RD model, while in Sec. III we briefly outline ε>0. Noticethatwhenanendreptonmovestowardsan
the basic ideas of DMRG. Sec. IV is dedicated to the emptylatticesite(tube renewalprocess)ithasintotald
properties of reptation in absence of an external field, in different possibilities of doing so forward and backward
particular to the scaling properties of τ and D. Sec. V inthefield(seeFig. 1),thereforetheassociatedratesare
collects a series of results of reptation in a field, while dB and dB−1 respectively. Onthe contrarythe end rep-
Sec. VI concludes our paper. ton can retract by moving to the (unique) site occupied
by its neighbor with rate B or B−1. Summarizing the
possible moves are: (a) stored length diffusion for inner
segments,i.e. (0, 1 1,0)withratesB andB−1,(b)
II. THE RUBINSTEIN - DUKE MODEL contractions for e±xter↔na±l reptons ( 1 0) with rates B
andB−1 and(c)stretchesforexter±nal→reptons(0 1),
with rates dB and dB−1. Notice that the par→am±eter
In the RD model the polymer is divided into N units, d enters in the model only at the stretching rates (c).
called reptons, which are placed at the sites of a d- Ratherthanlinkingdtothelatticecoordinationnumber
dimensionalhypercubic lattice (see Fig. 1). The number we interpret it as the ratio between stretching rates and
of reptons that each site can accommodate is unlimited rates associated to moves of inner reptons. This allows
and self-avoidance effects are neglected. Each configura- to choose any positive values for d.
tion is projected onto an axis along the (body) diagonal
of the unit cell and it is identified by the relative coor-
dinates y z z of neighboring reptons along the
i i+1 i
chain (zi in≡dicates−the projected coordinate of the i-th (cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)
repton). The relative coordinates can take three values
y = 1,0,1 and there are thus in total 3N−1 different
i
−
configurations for a chain with N reptons. (cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)
ε
8
5 (cid:0)(cid:1) (a) (cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1) (b)
(cid:0)(cid:1)
(cid:0)(cid:1)
(cid:0)(cid:1) 6 (cid:0)(cid:1)
4(cid:0)(cid:1) (cid:0)(cid:1)7 FIG.2. (a) In thick line we indicate a test chain moving
(cid:0)(cid:1) 3 (cid:0)(cid:0)(cid:1)(cid:1) inanenvironmentofotherchains. Bymodifyingitsends,for
(cid:0)(cid:1)
(cid:0)(cid:1) instance, by attaching molecules of large size one can lower
(cid:0)(cid:1) the parameter d, as in this way one expects that stretchings
1 (cid:0)(cid:1)2
(cid:0)(cid:0)(cid:1)(cid:1) of thechain out of thetubewill occur tolower rate. (b)The
(cid:0)(cid:1)
lowering of d would also occur for reptating polymers with
short branchingends.
FIG. 1. Example of a configuration in the RD model in
d = 2 dimensions. In terms of the relative coordinates
As we will see,a variationof d has some importantef-
projected along the field direction the configuration reads
fectsonthecorrectionstotheasymptoticbehavior. This
y={1,0,1,1,−1,1,0}. Anon-zerofield(ε)biasesthemotion
of the reptons, which occurs with rates B =exp(ε/2) (B−1) regime may turn out to be relevant for the experiments.
formovesinthedirection(oppositeto)theappliedfield. The Moreoverit is conceivablethatthe parameterdcouldbe
arrowsindicatethepossiblemoves. Noticethatanendrepton somewhatmodified in anexperiment. Figure 2 schemat-
canstretchtodlatticepositionsforwardandbackwardinthe ically illustrates a possibility of doing so, namely by at-
field (as repton 8 in theexample). taching to the chain ends some large molecules so that
2
for each chain tube renewal moves would be impeded. trivial (see next section). The right eigenvector gives
The same effects would be possible if the linear chain is the stationaryprobability distribution. The next to low-
modified suchas to havebranches close to the endpoints est eigenvalue, the gap 1/τ, yields the slowest relaxation
(Fig. 2(b)), or an end part stiffer than the rest of the time τ of the decaytowardsthe stationarystate. Gener-
chain (as it could be realized in block copolymers). In allythe DMRGmethodworksbestwhentheeigenvalues
both cases we expect that stretches out of the confining are well separated. For long chains and stronger driving
tube would be suppressed while chain retractions would fields the spectrum of H gets an accumulation of eigen-
not be much impeded. values near the zero eigenvalue of the stationary state.
Once the rates for elementary processes are given, the This hampers the convergence of the method seriously
stationary properties of the system can be found from and enlarging the basis m becomes of little help, while
the solution of the Master equation standardly this improves the accuracy substantially. In
order to construct the reduced density matrix from the
dP(y,t)
= Hyy′P(y′,t) (1) lowest eigenstates one needs to diagonalize the effective
dt − matrices H at each DMRG step, we used the so-called
y′
X
Arnoldimethodwhichis knowntobe particularlystable
inthelimitt . HereP(y,t)indicatestheprobability for non-hermitian problems [23].
→∞
of finding the polymer in a configurationy at time t and
the matrix H contains the transition rates per unit of
time between the different configurations of the chains, IV. ZERO FIELD PROPERTIES
as given in the rules discussed above.
Inthis sectionwepresenta seriesofDMRG resultson
the scaling behavior as function of the polymer length
III. DENSITY MATRIX RENORMALIZATION
N of the so-called tube renewal time τ(N), i.e. of the
characteristictime for reptation, and of the diffusion co-
DMRG was introduced in 1992 by S. White [18] as an
efficient D(N).
efficient algorithm to deal with a quantum hamiltonian
for one dimensional systems. It is an iterative basis–
truncation method, which allows to approximate eigen-
A. Tube renewal time
values and eigenstates, using an optimal basis of size m.
It is not restricted to quantum system; it has also been
successfullyappliedtoaseriesofproblems,rangingfrom The renewal time, τ, is the typical time of the repta-
two dimensional classical systems [21] to stochastic pro- tion process, i.e. the time necessary to loose memory of
cesses [22]. The basic common feature of these problems any initial configuration through reptation dynamics. τ
is the formal analogy with the Schr¨odinger equation for is given by the inverse of the smallest gap of the matrix
a one dimensional many–body system. DMRG can also H,which canbe seen asfollows: Starting fromaninitial
be applied to the Master Equation (1) for the reptat- configurationonehas,asymptoticallyforsufficientlylong
ing chain, albeit with some limitations due to the non– times (t ):
→∞
hermitiancharacterofthematrixH. Inthisspecificcase
westartfromsmallchains,forwhichH canbediagonal- P(t) = φ0 +e−t/τ φ1 +... (2)
| i | i | i
ized completely, and construct effective matrices repre-
senting H for longer chains through truncation of the where φ0 and φ1 are eigenstates of the matrix H de-
| i | i
configurational space followed by an enlargement of the fined in Eq. (1) with eigenvalues E0 = 0 and E1 = 1/τ,
chain. Thisisdonethroughtheconstructionofareduced respectively. The state φ0 is the stationarystate of the
| i
densitymatrixwhoseeigenstatesprovidetheoptimalba- process, while the gap of H corresponds to the inverse
sis set, as can be shown rigorously (see Refs. [18,19] for relaxation time. Notice that, as the matrix H is sym-
details). The size m of the basis remains fixed in this metric at zero field, E1 is always a real number, while
process. By enlarging m one checks the convergence of for non-Hermitian matrices the eigenvalues may get an
the procedure to the desired accuracy. Hence m is the imaginarypart, whichcausesa relaxationto equilibrium
main control parameter of the method. In the present through damped oscillations.
case we found that m=27 is sufficient for small driving An important point of which we took advantage in
fields and we kept up to m=81 state for strongerfields. the calculation is the fact that the ground state φ0 of
| i
H is only hermitian for zero driving field. So for a fi- the matrix H is known exactly in absence of any driv-
nite driving field one needs to apply the non–hermitian ing fields. Since H is stochastic and symmetric (i.e.
variantofthestandardDMRGalgorithm[22]. Foranon– y′Hyy′ = y′Hy′y = 0), it follows immediately that
hermitian matrix H one has to distinguish between the the stationary state is of the form φ (y) = c, with c a
0
P P
right and left eigenvector belonging to the same eigen- constant. Instead of applying the DMRG technique di-
′
value. Since H is a stochastic matrix the lowest eigen- rectly to the matrix H we applied it to the matrix H
value equals 0 and the corresponding left eigenvector is defined by:
3
′
H =H +∆φ φ . (3) uncertainties,as these are amplifiedwhen taking numer-
0 0
| ih |
ical derivatives. We stress that already in the log-log
′
Now if we choose ∆> E1 the lowest eigenvalue of H is plot of Fig. 3 the deviation of the data from linearity
equaltoE1,thereforethe problemofcalculatingthe gap is noticeable, therefore the values of z given above are
of H is reduced to the calculation of the ground state of justaveragevaluesandnottobeexpectedasasymptotic
′
a new (non-stochastic)matrixH . This approachis con- ones.
siderably more advantageous in term of CPU time and Figure 4 shows a plot of z for the data given in Fig.
N
memory required for the program (for more details see 3, plotted as function of 1/√N. In the thermodynamic
Ref. [24]). limit N , z is seen to converge towards z = 3,
N
Figure 3 shows a plot of lnτ(N) vs. lnN for the re- in agreem→ent∞with de Gennes’ theory [1]. Corrections to
newal time as calculated from DMRG methods for vari- this asymptotic limit yieldz >3 whenN is sufficiently
N
ous values of the stretching rate d and for lengths up to large, with deviations towards z < 3 for small d and
N
N =150. Thedatahavebeenshiftedalongtheabscissae not too large N.
axis by an arbitrary constant. We fitted the data by a
linear interpolation, as it is mostly done in Monte Carlo 4
simulationsandinexperiments. Theresultingslopespro-
CLF
vide estimates of the exponent z, which we find to be
sensitive to the stretching rate d. For d sufficiently large
d=3
(d > 2) we find z 3.3, in agreement with experiments
[4] and previous M≈onte Carlo simulation results [11,25]. 3.5
The value of z decreases when d is decreased. For d suf-
d=1
ficiently small (d 0.1) we find z 2.7 < 3, which is a N
≈ ≈ z
regime that has not yet been observed in experiments.
3
20
d=0.25
d = 0.1, z = 2.70
d = 0.5, z = 3.00
d = 1.0, z = 3.21 d=0.1
15 d = 3.0, z = 3.29
2.5
st. 0 0.1 0.2 0.3 0.4
on N−1/2
c
+ 10
)
N FIG. 4. Solid lines: Effective exponent zN for various d.
(
τn Dashed line: Corrections due to CLF as predicted by Doi
l
(Eqs. 5 and 6) for the d = 3 curve. Dot-dashed line: In-
5
clusion of higher order corrections as given by Eq. 8 for the
curved=1.
0 Thereexistsometheoreticalpredictionsfortheformof
1 2 3 4 5 thecorrectionstotheasymptoticbehavioroftherenewal
ln N time[26]. Thesearebasedoncontourlengthfluctuations
(CLF), i.e. on the idea that the length of the tube fluc-
FIG.3. Plot of lnτ(N) vs. lnN for various d. We report tuates in time and that this would help accelerating the
the values of the slope of the data obtained from a linear renewal process. Such fluctuations were not taken into
interpolation. account in the originalwork of de Gennes, who assumed
the tube to have a fixed length. According to CLF the-
To shed some more light into these results we consid- ory, thus, τ scales as [26]
ered the effective exponent, i.e. the following quantity:
2
N
lnτ(N +1) lnτ(N 1) τ(N) N3 1 0 , (5)
zN = ln(N +1)−ln(N −1) , (4) ∼ −rN !
− −
withN atypicallength. Substituting this equationinto
0
which is the discrete derivative of the data in the log-log
Eq. (4) one obtains:
scale plot. Such quantity probes the local slope (at a
givenlengthN)ofthedataofFig. 3andprovidesabet- N /N
(CLF) 0
z =3+ . (6)
terestimateofthefiniteN correctionstotheasymptotic N
1 N /N
p 0
behavior. The calculation of z requires very accurate −
N
dataandcannotbeeasilyperformedbyMonteCarlosim- which is a monotonically increpasing function of 1/√N.
ulation results which are typically affected by numerical For all lengths in the physical range N > N one has
0
4
z(CLF) > 3. Eq. (6) is plotted as a dashed line in Fig. B. Diffusion coefficient
N
4, where we have chosen N = 2.3 in order to obtain
0
the best fit of the DMRG data with d = 3. This value
is consistent with that predicted by the CLF theory [3].
Weconsidernowasimilarscalinganalysisofthediffu-
OurresultsfurthersuggestthatthevalueofN increases
0 sion coefficient D. To compute the diffusion constant as
when d is decreased.
functionofthe chainlengthweusedtheNernst–Einstein
While for sufficiently long chains the data apparently
relation [14]:
merge to the CLF theory given by Eq.(6) higher order
corrections appear to be of opposite sign. When d is
v
loweredthelatterbecome particularlystrongsothatthe D = lim . (9)
ε→0Nε
effective exponent for a certain range of lengths is even
smaller than 3. To our knowledge such an effect has
not yetbeen observedin experiments. Presumablystan- where v is the drift velocity of the polymer with N rep-
dard polymer mixtures will correspond to d > 1, which tons and subject to an external field ε. In order to esti-
is a regime where an exponent z 3.4 is ob∼served [4]. mateD(N)weconsideredsmallfields(downtoε 10−5)
It is conceivable that mixtures wi≈th modified architec- and calculated the velocity vi of the i-th repton∼in the
ture, as those illustrated in Fig. 2(b), i.e. long polymers stationarystate using the formulas givenin SectionVB.
with short branching ends, would correspond to curves As a check for the accuracy of the procedure we verified
at much lower d. Experimental results for such systems that the drift velocity vi is independent on i, i.e. the
wouldbeofhighinterestinordertocheckthepredictions positionofthe reptonalongthechaininwhichitis com-
of the RD model in the low d regime. puted.
In order to gain some more insight in the precise form
TABLEI. Comparison between diffusion coefficients from
of τ(N) we considered higher order terms which lead to
exact diagonalization methods, as given in Ref. [17], and as
an expression of the type:
obtainedfromthecalculationofr(ε)=v/NεfromtheDMRG
algorithm with m=18 for decreasing ε.
2
N N
τ(N) N3 1 0 +A 0 (7)
∼  −rN ! N  N r(ε=10−2) r(ε=10−3) r(ε=10−4) DN2 (Ref. [17])
5 0.864873 0.864907 0.864908 0.864908
 
7 0.769024 0.769057 0.769057 0.769057
Recently, Milner and McLeish [27] formulated a more
9 0.703975 0.703951 0.703951 0.703951
complete theorybeyond that ofCLF by Doi, using ideas
11 0.657701 0.657550 0.657549 0.657549
fromthetheoryofstressrelaxationforstarpolymers[28], 13 0.623372 0.623010 0.623006 0.623006
which to lowest orders in 1/√N yields an expression of 15 0.596975 0.596304 0.596297 0.596297
the type given above. The effective exponent now reads: 17 0.576097 0.575007 0.574996 0.574996
19 0.559224 0.557595 0.557580 0.557579
N 1 (1+A) N /N
0 0
z =3+ − (8)
N N 2
r 1 N /N p+AN /N
0 0
−
5
(cid:16) p (cid:17)
Notice that the CLF expression (Eq. 6) diverges for d = 0.1, x = 1.80
N N , while the previous formula this divergence d = 0.5, x = 1.80
0
→
does not occur. For N not too large the numerator in d = 1.0, x = 2.10
Eq. (8) may change sign, reproducing features found in d = 3.0, x = 2.12
the DMRGcalculations,i.e. aneffectiveexponentz <3. st. 0
n
Thedot-dashedlineofFig. 4representsafitforthed=1 o
c
case using Eq.(8); we find that the choice N0 = 3.6 and N)+
A = 0.44 fits very well the numerical data for N > 25, (
D
whiledeviationsforshorterchainsareclearlyvisible. We n −5
l
stress that Eq. (8) fits quite well the renewal time data
in the cage model [29], which is another lattice model
of reptation dynamics [31]. Finally, very recently the ef-
fect of constraints release (CR), introduced in the RD
model in a self-consistent manner, has been considered −10
[30]. Apartfromsmallquantitativeshifts inthe effective 1 2 3 4 5
ln N
exponents, the is unchanged.
FIG. 5. Log-log plot of the diffusion constant as function
ofthechainlengthforvariousvaluesofthestretchingrated.
5
In practice we estimated the ratio r = vN/ε for de- recent experimental investigation of diffusion coefficient
creasing ε until convergence was reached (typically we for polymer melts and solutions yielded some different
considered ε 10−4 10−5). Table I shows the numeri- results concerning its scaling behavior. Early measure-
calvaluesofr≈(ε)forv−ariousεobtainedfromDMRGwith ments of polymer melts yielded D N−2 [5], while in
m=18 states kept, for short chains,so that they canbe concentrated solutions typically D ∼N−2.4 [6,7]. Very
∼
compared with the exact diagonalization data shown in recently, however, on hydrogenated polybutadiene con-
the last column. The results for ε = 10−4 are in excel- centrated solutions and melts it was found D N−2.4
lent agreementwith exact diagonalizationdata for DN2 for both cases [9]. A reanalysis of previous expe∼rimental
(from Ref. [17]) reported in the last column, which indi- results on several different polymers lead to the conclu-
catesthattheprocedureusedtoextrapolatethediffusion sion D N−2.3 [8]. Thus the issue experimentally has
∼
coefficient from the Nernst-Einstein relation (Eq. 9) is not yet fully settled. A theoretical analysis of the con-
correct and reliable. tour length fluctuations on the diffusion coefficient was
According to reptation theory the diffusion coefficient recently performed [34] leading to the expression:
for large N should scale as D(N) 1/N2 [1]. This re-
∼
sults has been also derived rigorously for the RD model −1
where the coefficient of the leading term is also known D (N) N−2 1 N0 (13)
[32,15]: CLF ∼ −rN !
1
D(N)= (10)
(2d+1)N2 In Fig. 6 we plotted the corresponding effective expo-
nent, fitted to the d = 3 results. Again, the numerical
Early experimental results on diffusion coefficients for
results seem to approach the CLF formula only for very
polymermeltswereconsistentwithapower 2[5],while
− longchains,wheretheeffectiveexponentisalreadyquite
morerecentresultssuggestthat suchexponentwouldbe
close to the asymptotic value. Our results indicate some
significantlyhigher[8]. Infacttheoriginalmeasurements
variation of the effective exponent as function of the pa-
consistent with a scaling 1/N2 raised quite some prob-
rameterd. Itwouldbethereforeinterestingtoinvestigate
lems in the past, which was pointing to an inconsistence
the scaling of D for the polymer with short branching
between the scaling of D and τ due to the following ar-
ends,asthoseillustratedinFig. 2(b)inordertotestthe
gument. For a reptating polymer, D and τ the typical
validity of our predictions at smaller d.
radius of gyration should scale as:
2.5
R √Dτ (11)
g ∼ CLF
2.4
Now as the polymer obeys gaussianstatistics R √N, d = 3.0
g ∼ 2.3
which implies that, if τ N3.4 then necessarily D
N−2.4, for the same range∼of lengths. This argument∼is 2.2
however not quite correct in this form since R √N d = 1.0
g 2.1
∼
only asymptotically, thus finite size corrections may af-
)
fect D and τ differently, as indeed happens in the RD (N 2
x
model. 1.9 d = 0.5
Figure 5 shows a plot of lnD as a function of lnN for
1.8
various values of the stretching rate d. The best fitting
parametersx for the scalingD 1/Nx aregiven. As for 1.7
∼
the diffusion constantthe exponent passes the presumed d = 0.1
1.6
asymptotic value of 2 for d = 1 and d = 3, while x < 2
for d = 0.5 and d = 0.1. Other numerical investigations 1.5
0 0.1 0.2 0.3 0.4 0.5
oftheRDandrelatedmodelyieldedx 2.0[11],x 2.5
≈ ≈ N−1/2
[25], and x 2.0 [33]. Again, it is best to analyze the
≈
effective exponent
FIG.6. Effective exponent xN as defined by Eq. (12) for
various d. The dotted line is the theoretical prediction from
lnD(N +1) lnD(N 1)
x = − − , (12) Eq. (13).
N
− ln(N +1) ln(N 1)
− −
which is shown in Fig. 6. x shows a similar behav- Finally we mention that the further order correction
N
ior as the renewal exponent. However, comparing the terms for the scaling of D have raised recently some
same values of N and d in Figs. 4 and 6 one notices debate about the nature of the expansion [16,15]. The
that finite N corrections are weaker for the diffusion co- DMRGresults,whichareaccurateenoughto investigate
efficient. For instance, for d = 3 and N−1/2 = 0.3 one the higher orders, clearly show that these results can be
has z 3.8, while x 2.3. As mentioned above, very well represented by a series in 1/√N [20].
N N
≈ ≈
6
V. REPTATION IN THE PRESENCE OF AN 0.14
ELECTRIC FIELD 0
0.12 v
n
l −10
Forfinitevaluesofǫseverelimitationsappearprevent-
0.10
ingustoextendthecalculationstothelongchainswhich
we were able analyze in the limit ǫ 0. These limita- 0 5
→ 0.08 ε
tions are intrinsic to the problem and the Arnoldi rou-
tine for finding the lowest eigenvalue of the matrix H. v
For finite ǫ and longer chains a massive accumulation of 0.06
smalleigenvaluesnearthe groundstate eigenvaluestarts
to emerge. The problem is also present in the straight 0.04 N=5
N=7
application of the routine to small chains, for which an
N=9
exactdiagonalizationofthe matrixH canbe performed. 0.02 N=11
For N >11 and ǫ 1 the method for finding the lowest N=15
≥
eigenvalue is no longer convergent. The DMRG method, 0.00
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
describedabove,yieldsconvergentresultsuptothechain ε
N = 15. For smaller fields, ε < 1, which is however the
most interesting region from the physical viewpoint, the
FIG.7. The velocities as a function of an electric field for
calculation can be extended to somewhat longer chains
chains of different lengths for d = 1. Inset: lnv vs. ε for
i.e. N 30 40. Neither extension of the truncated
≈ − N =7;theasymptoticbehaviorisconsistentwiththepredic-
basis set, nor the inclusion of more target states is a
tion from Eq.(14), shown as a dashed line.
remedy for the lack ofconvergence. Despite the intrinsic
difficulty oftreating long chains,substantialinformation
We have also investigated the behavior of the velocity
can be obtained by the analysis of the drift velocity for
as function of the stretching ratio d. In Fig. 8 we have
intermediate chain lengths as function of ε.
plotted the various curves for chain length N = 7 for d
rangingfrom 0.05 to 3.25. We see that the overallveloc-
ity becomes small for small d which is simply a sign of
A. Drift velocity vs. field
slowing down due to the slowed down motion of the end
reptons. Forthehighervaluesofdweseeagainadecrease
In Fig. 7 we have plotted the drift velocity as func- in the drift velocity, which is explained by stretching of
tion of ǫ for chains up to N = 15 and stretching ratio the chain as we will see.
d = 1. The general behavior is a rise turning over into
an exponential decay. For small ε the rise is linear, in 0.20
agreement with the results presented for the zero field d = 0.05
diffusioncoefficient. Infactwe havecalculatedthe diffu- d = 0.50
d = 1.16
sioncoefficientasthederivativeofthedriftvelocitywith
0.15
d = 1.50
respect to ε (see Sec. (IVB)). The exponential decay is
d = 2.00
inaccordancewiththebehaviorpredictedbyKolomeisky
d = 3.25
[35]. For largefields the probabilitydistributionnarrows
down to a single U shaped configuration with equally v 0.10
long arms at both sides of the center of the chain. For
odd chains Kolomeisky shows that, in the limit ε
→ ∞
[35] 0.05
2 N
v exp − ε (14)
∼ 2
(cid:18) (cid:19) 0.00
0 1 2 3 4 5
ε
For the chains up to N = 7 we could confirm this be-
havior(see inset ofFig. 7), but forthe longerchains one
wouldhavetogotolargervaluesofǫthancanbehandled FIG.8. The velocities as a function of an electric field for
withtheArnoldimethodtoseetheasymptoticbehavior. different d. The numberof reptons is fixed at N =7.
Generally the curves area compromisebetween twoten-
dencies: the expression for the drift grows with the field Acloserinspectionofthe curvesv vs. εinthe interval
buttheprobabilityonaconfigurationshiftstowardscon- 0<ε 1forthelongerchainsshowsadefinitedeviation
≤
figurations with the leastvelocity. With increasing ǫ ini- ofthelinearrisetolargervaluesofthedriftvelocity. This
tially the first effect dominates but gradually the second interesting intermediate regime was also investigated by
effect overrules the first. Barkemaetal.[16]bymeansofMonteCarlosimulations.
7
They propose to describe the drift velocityin this region the context of the biased reptation model (BRM) [37],
by the phenomenological crossoverexpression whichprovidesasimplifiedpictureofthephysicsofpoly-
ǫ mer reptation in a field, in which fluctuations effects are
v(ǫ,N) [1+A(ǫN)2)]21 (15) neglected. More refined analytical calculations in mod-
≃ (2d+1)N
els where fluctuation effects are taken into account still
which seems to fit quite well the Monte Carlo data [36]. predictthepresenceofavelocityminimum[38],thusthe
Thus in the regime where ǫ is small but the combination original result of the BRM was confirmed.
ǫN oforderunity,thedriftvelocitybecomesindependent It should be pointed out that the Eq. (15) does not
of N (band collapse) and quadratically dependent on ǫ. yield any minimum in the velocity as ∂v/∂N < 0 for
all N, thus it cannot be strictly correct. Although our
0.060 results are limited to the onset of the minimum, we sus-
0.0105 ε=0.01 pect that the minimum is not a very pronounced one so
ε=0.05
0.050 ε=0.10 thatforlongerchainsthecurvesv vs. 1/N wouldappear
ε=0.15 rather flat. A shallow minimum could have been thus
0.0100 ε=0.20 missed in the Monte Carlo simulations of Ref. [16].
0.040 ε=0.25
ε=0.30 Indeed, the data presented in Ref. [16] yield, for in-
0.0095 stance, v 0.010 for ε = 0.2 and N = 100, thus the
v 0.030 0.03 0.04 0.05 0.06 0.07 limiting ve≈locity appears to be approached very slowly
from below (see inset of Fig. 9). Other previous simu-
0.020
lations by Duke [12] of the RD model yielded instead a
quite clear velocity minimum, but it is difficult to judge
0.010 the quality of his data as no error bars are reported.
Probably the inversion effect is somewhat weaker than
0.000 reported in [12]. We investigated further the effect of
0.00 0.05 0.10 0.15 0.20 d and found that for higher d the minimum appears to
1/N
be somewhat more pronounced, which may explain the
results of Ref. [12], where d = 6 was taken (the lattice
FIG. 9. Plot of the drift velocity v as function of the in- coordination number for an fcc lattice).
versechainlength1/N,forvariousfixedvaluesoftheelectric Unfortunately, present DMRG computations are lim-
field ε and d=1. Inset: Blow up of v vs. 1/N for ε=0.2; a itedtotheonsetoftheinversionphenomenon. Theprob-
minimum in thevelocity can beclearly seen. lem with the DMRG approach is that it builds up an
optimal basis for the stationary state using information
In Figure 9 we plotted the velocity as function of the of stationary states for shorter chains. Around the band
inverse chain length 1/N for various values of the elec- inversion point there is a kind of phase transition from
tric field. The deviations from the linear regime, where short unoriented polymers to long oriented ones (see for
v ε/N, are clearly observable as for sufficiently long instance Ref. [10]). This change implies that the opti-
∼
chains and not too small fields the drift velocity ap- mal basis for short chains may be no longer a good one
proaches a non-vanishing value. A closer inspection on when longer chains are considered. We alleviated some-
the curves reveals that the limiting constant velocity in what this problem using the DMRG algorithm at fixed
theasymptoticregimeN isreachedthroughamin- Nε, which allows indeed to study slightly longer chains,
→∞
imum in the velocity ata givenpolymer length Nmin. In yetnotenoughtogodeepinto the bandcollapseregime,
our calculations this occurs for any values of the fields butsufficienttobringinevidencethe velocityminimum.
in the range ε 0.10. Presumably the minimum is
Originally,itwasthoughtthatthedriftvelocityshould
≥
shifted to lengths N much longer than those which can
be described by a scaling function in terms of the com-
be reachedin the presentcalculation. The existence of a bination Nε2 [39,40]. It is nowadays believed that the
minimum velocity at finite length can be clearly seen in
correct scaling form for the velocity should be given by
the insetofFig. 9whichplotsv vs. 1/N forε=0.2,but
an expression [41,16]
it is present also in other curves. According to our er-
rorestimatestypicaluncertaintiesinthevelocitiesareat ε
v(ε,N)= g(Nε) (16)
thesixthdecimalplace,thuserrorbarsaremuchsmaller N
than all symbol sizes in the figure and inset.
Thepossibleexistenceofavelocityminimumhasgiven with g(x) g > 0 for x 0 and g(x) x for x
0
→ → ∼ → ∞
rise to quite some discussionsin the pastandit is a phe- in order to match the known behavior of the velocity at
nomenon referred to as band inversion (for a review see largeandsmallfields. The conditionto havebandinver-
Ref. [10]). In the band inversion region long polymers sion is ∂v/∂N = 0 for some N at fixed fields, which is
′
migrate faster along the field direction than short ones, equivalent to the requirement g(x) = xg (x), for a non-
a physical situation that looks at first sight somewhat zerovalueofx=Nε. Noticethatchoosingthemostgen-
counterintuitive. Band inversion was predicted first in eral scaling form v(ε,N) = ε g(Nεα) and from the re-
N
8
′
quirement of ∂v/∂N =0 one still obtains g(x)=xg (x), p (y)=1. (19)
i
with x=Nεα. y
X
The scaling behavior of the minimum as function of
the field may be used to test the velocity formula. We
′
The9possiblevaluesofp (y,y )arerestrictedbythese
calculated the polymer length N at which the mini- i
min
conditions. One finds another set of relations between
mumofthevelocityoccursasfunctionoftheappliedfield
these quantities by summing the Master Equation over
ε. If Eq. (16) is correct we expect ε 1/N . Figure
min
∼ all but one segment value. Exclusion of an internal seg-
10 shows a plot of lnε vs. lnN for d = 1. The data
min
ment y from the summation yields for y =y = 1
show some curvature due to corrections to scaling and i i ±
approach,forthelongestchainsanalyzed,theasymptote
ε N−α, with α = 1.4, as illustrated in the figure. As Bypi−1(0,y) B−ypi−1(y,0)=
Nm∼in is increased one observes a systematic decrease of Bypi(0,y) B−−ypi(y,0)=v(y). (20)
thelocalslopeofthedatasuggestingthattheasymptotic −
exponentα<1.4. Inanattempttoreachtheasymptotic
Note that v(y) is independent of the index i of the seg-
regimeweextrapolatedthelocalslopesofthedatainFig.
ment under consideration. The two relations (20) are an
10 in the limit N , which yield an extrapolated
min
→ ∞ expressionof the fact the averagevelocity of the reptons
value α 1.1, not far from the prediction of Eq.(16),
≈ inthefielddirectionandthecurvilinearvelocityarecon-
but not fully consistent with it. To prove Eq.(16) more
stant along the chain. Taking the segment value y = 0
convincinglyonewouldneedtoinvestigatelongerchains. j
one finds
−1.0
v =v(1) v( 1)=2v(1), (21)
− −
DMRG
Slope=1 with v the drift velocity. Excluding the end segments
−1.5 Slope=1.4 from the summation yields the 2 equations (y = 1)
±
εn v(y)=Byp1(y) dB−yp1(0)
l v(y)=dByp (0)− B−yp (y) (22)
N N
−
−2.0
One observes that the 3 probabilities on the end seg-
mentsarefixedbythenormalizationandthe2equations
(22). Ingeneraltheserelationsshowhowdelicatethede-
−2.5
velopment of the correlations is. In the fieldless case the
−3.6 −3.4 −3.2 −3.0 −2.8
ln(1/N ) probabilityofaconfigurationfactorizesinaproductover
min
probabilities of segments. So for d = 1 the probability
on any segment becomes equal to 1/3. Then of course
FIG.10. Plot of lnε vs. lnNmin for d=1. the velocitiesvanish. Consideringthe termslinearinthe
field (or in B B−1) one sees that the drift velocity of
−
the first repton, given by
B. Correlations and profiles
v =Bp (1) B−1p ( 1)+d(B B−1)p (0), (23)
1 1 1
− − −
A more detailed insight in the shapes of the config-
urations is obtained by plotting averages of the local
requires a delicate compensation in the linear deviations
variables y . The DMRG procedure leads naturally to
i of the probabilities p (y) in order to give a value which
1
the determination of the probabilities for two consecu-
vanishes as 1/N for long chains.
tive segments
Rather than giving the values of p (y) we plot the av-
i
′ erages
pi(y,y )=hδy,yiδy′,yi+1i. (17)
The probabilities on a single segment follow from these yi =pi(1) pi( 1) (24)
h i − −
values
and
′ ′
pi(y)= pi(y,y )=hδy,yii= pi−1(y ,y), (18)
Xy′ Xy′ yiyi+1 =pi(1,1)+pi( 1, 1) pi(1, 1) pi(1 ,1).
h i − − − − − −
which in turn are normalized (25)
9
5 6
2
5
4
> >
410<y>l 0 <yykk+1 10<y>l 0 <yykk+1
dd==00..1205 610 εε==00..00051 102
−5 d=0.50 ε=0.24
d=1 ε=0.50
d=2 ε=0.70
d=3
0 −5 0
−1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0
l (i) k (i) l (i) k (i)
N N N N
FIG. 11. The profiles for ε = 0.001 plotted as function FIG. 13. As in Fig. 11 for d = 1 and varying ε. The
of the reduced distances lN(i) = (2i − N)/(N − 2) and numberof reptons is fixedat N =15.
kN(i)=(2i−N +1)/(N −3) (with i=1,2,...N −1). The
numberof reptons is fixed at N =15.
As an example of the behavior in the intermediate
regime(ofthevelocityprofilesFig. 3)wehaveplottedin
Fig. 13 the situation for N =15, d=1 and various val-
ues of ǫ. For the larger values of ǫ the average y does
0 h ii
notchangeverymuchbuttheplateauofthecorrelations
y y in the middle keeps rising. We expect that for
i i+1
h i
longerchainsalargeregioninthemiddledevelopswhere
−10
theshapevariesweaklywithǫandwherethecorrelations
between consecutive segments increase. Thus the chain
>
yi obtains longer stretches which are oriented in the field,
<
<j −20 either up or down, but which largely compensate, such
Σi
4 0 d=0.10 thattheoverallshapeinthe middleremainsmoreorless
1 d=0.25 the same. This is in agreement with the speculations of
−30 d=0.50 Barkema et al. [16] on the chain as a stretched sequence
d=1.00 of more or less isotropic ”blobs”.
d=2.00
d=3.00
−40 10
0 5 10 15 ε=0.05
j ε=0.50
ε=1.16 8
6 ε=1.50
ε=2.00
FIG.12. Theshapesofthechainforε=0.001andvarious
ε=2.50
values of d. The numberof reptons is fixed at N =15. ε=3.25
2
ε=4.00 >
<y>l yk+1
A typical plot for the small ǫ regime is given in Fig. 0 yk
1 <
11. The globalsymmetry due to the interchange of head −2 10
and tail makes the average(24) anti–symmetric with re-
specttothemiddleandtheaverage(25)symmetric. One
−6
observesthatthe featuresincreasewiththemobilitydof 0
the end reptons. The averages give information about
the average shape of the chain. In Fig. 12 we translate
−10
the averages (24) into average spatial configurations by −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0
integrating (summing) the segment values to positions l (i) k (i)
N N
with respect to the middle repton. Clearly the develop-
ment of the U shape is visible with increasing d. The FIG.14. The profiles for d=1. The number of reptons is
effect is larger at the ends than in the middle. fixedat N =9.
10