Intramolecular long-range correlations in polymer melts:
The segmental size distribution and its moments
J.P. Wittmer,∗ P. Beckrich, H. Meyer, A. Cavallo, A. Johner, and J. Baschnagel†
Institut Charles Sadron, CNRS, 23 rue du Loess, 67037 Strasbourg Cédex, France
(Dated: November 1, 2018)
Abstract
Presenting theoretical arguments and numerical results we demonstrate long-range intrachain
correlations in concentrated solutions and melts of long flexible polymers which cause a systematic
swelling of short chain segments. They can be traced back to the incompressibility of the melt
leading to an effective repulsion u(s) ≈ s/ρR3(s) ≈ ce/
s when connecting two segments together
where s denotes the curvilinear length of a segment, R(s) its typical size, ce ≈ 1/ρb3e the “swelling
coefficient”, be the effective bond length and ρ the monomer density. The relative deviation of
the segmental size distribution from the ideal Gaussian chain behavior is found to be proportional
to u(s). The analysis of different moments of this distribution allows for a precise determination
of the effective bond length be and the swelling coefficient ce of asymptotically long chains. At
striking variance to the short-range decay suggested by Flory’s ideality hypothesis the bond-bond
correlation function of two bonds separated by s monomers along the chain is found to decay
algebraically as 1/s3/2. Effects of finite chain length are considered briefly.
PACS numbers: 61.25.Hq,64.60.Ak,05.40.Fb
∗Electronic address: jwittmer@ics.u-strasbg.fr
†URL: http://www-ics.u-strasbg.fr/~etsp/welcome.php
http://arxiv.org/abs/0704.1620v2
mailto:jwittmer@ics.u-strasbg.fr
http://www-ics.u-strasbg.fr/~etsp/welcome.php
I. FLORY’S IDEALITY HYPOTHESIS REVISITED
A cornerstone of polymer physics. Polymer melts are dense disordered systems con-
sisting of macromolecular chains [1]. Theories that predict properties of chains in a melt
or concentrated solutions generally start from the “Flory ideality hypothesis” formulated
already in the 1940s by Flory [2, 3, 4]. This cornerstone of polymer physics states that
chain conformations correspond to “ideal” random walks on length scales much larger than
the monomer diameter [1, 4, 5, 6]. The commonly accepted justification of this mean-field
result is that intrachain and interchain excluded volume forces compensate each other if
many chains strongly overlap which is the case for three-dimensional melts [5]. Since these
systems are essentially incompressible, density fluctuations are known to be small. Hence,
all correlations are supposed to be short-ranged as has been systematically discussed first
by Edwards who developed the essential statistical mechanical tools [6, 7, 8, 9, 10] also used
in this paper.
One immediate consequence of Flory’s hypothesis is that the mean-squared size of chain
segments of curvilinear length s = m − n (with 1 ≤ n < m < N) should scale as R2e(s) ≡
〈r2〉 = b2es if the two monomers n and m on the same chain are sufficiently separated along
the chain backbone, and local correlations may be neglected (1 ≪ s). For the total chain
(s = N−1 ≫ 1) this implies obviously that R2e(N−1) = b2e(N−1) ≈ b2eN . Here, N denotes
the number of monomers per chain, r the end-to-end vector of the segment, r = ||r|| its
length and be the “effective bond length” of asymptotically long chains [6]. (See Fig. 1 for
an illustration of some notations used in this paper.) For the 2p-th moment (p = 0, 1, 2, . . .)
of the segmental size distribution G(r, s) in three dimensions one may write more generally
Kp(s) ≡ 1−
(2p+ 1)!
〈r2p〉
(b2es)
= 0 (1)
which is, obviously, consistent with a Gaussian segmental size distribution
G0(r, s) =
2πsb2e
. (2)
Both equations are expected to hold as long as the moment is not too high for a given
segment length and the finite-extensibility of the polymer strand remains irrelevant [6].
Deviations caused by the segmental correlation hole effect. Recently, Flory’s hypothe-
sis has been challenged both theoretically [11, 12, 13, 14, 15] and numerically for three-
dimensional solutions [16, 17, 18, 19, 20] and ultrathin films [21, 22]. These studies suggest
that intra- and interchain excluded volume forces do not fully compensate each other on
intermediate length scales, leading to long-range intrachain correlations. The general phys-
ical idea behind these correlations is related to the “segmental correlation hole” of a typical
chain segment [19]. As sketched in Fig. 2, this induces an effective repulsive interaction when
bringing two segments together, and swells (to some extent) the chains causing, hence, a
systematic violation of Eq. (1). Elaborating and clarifying various points already presented
briefly elsewhere [18, 19, 20], we focus here on melts of long and flexible polymers. Us-
ing two well-studied coarse-grained polymer models [23] various intrachain properties are
investigated numerically as functions of s and compared with predictions from first-order
perturbation theory. (For a discussion of intrachain correlations in reciprocal space see
Refs. [14, 15, 19].)
Central results tested in this study. The key claim verified here concerns the deviation
δG(r, s) = G(r, s) − G0(r, s) of the segmental size distribution G(r, s) from Gaussianity,
Eq. (2), for asymptotically long chains (N → ∞) in the “scale-free regime” (1 ≪ s ≪ N).
We show that the relative deviation divided by ce/
s scales as a function f(n) of n = r/be
δG(r, s)/G0(r, s)
= f(n) =
+ 9n−
. (3)
As we shall see, this scaling holds indeed for sufficiently large segment size r and curvilinear
length s. The indicated “swelling coefficient” ce has been predicted analytically,
24/π3
(ρ being the monomer number density), where we shall argue that the bond length of the
Gaussian reference chain of the perturbation calculation must be renormalized to the effective
bond length be. Accepting Eq. (3) the swelling of the segment size is readily obtained by
computing 〈r2p〉 = 4π
dr r2+2pG(r, s). For the 2p-th moment this yields
Kp(s) =
3(2pp!p)2
2(2p+ 1)!
. (5)
For instance, for the second moment (p = 1) this reduces to K1(s) = 1 − Re(s)2/b2es =
s ≈ ce/
s. We have replaced in Eq. (5) the theoretically expected swelling coefficient
ce by empirically determined coefficients cp. It will be shown, however, that cp/ce is close
to unity for all moments. Effectively, this reduces Eq. (5) to an efficient one-parameter
extrapolation formula for the effective bond length be of asymptotically long chains albeit
empirical and theoretical swelling coefficients may slightly differ. While we show how be may
be fitted, no attempt is made to predict it from the operational model parameters and other
measured properties such as the microscopic structure or the bulk compression modulus
[6, 9, 10].
Outline. We begin our discussion by sketching the central theoretical ideas in Sec. II.
There we will give a simple scaling argument and outline very briefly some elements of the
standard perturbation calculations we have performed to derive them (Sec. II B). Details of
the analytical treatment are relegated to Appendix A. The numerical models and algorithms
allowing the computation of dense melts containing the large chain lengths needed for a clear-
cut test are presented in Sec. III. Our computational results are given in Sec. IV. While
focusing on long chains in dense melts, we explain also briefly effects of finite chain size. The
general background of this work and possible consequences for other problems of polymer
science are discussed in the final Sec. V.
II. PHYSICAL IDEA AND SKETCH OF THE PERTURBATION CALCULATION
A. Scaling arguments
Incompressibility and correlation of composition fluctuations. Polymer melts are essen-
tially incompressible on length scales large compared to the monomer diameter, and the
density ρ of all monomers does not fluctuate. On the other hand, composition fluctuations
of labeled chains or subchains may certainly occur, however, subject to the total density
constraint. Composition fluctuations are therefore coupled and segments feel an entropic
penalty when their distance becomes comparable to their size [12, 19]. As sketched in
Fig. 2(a), we consider two independent test chains of length s in a melt of very long chains
(N → ∞). If s is sufficiently large, their typical size, R(s) ≈ be
s, is set by the effective
bond length be of the surrounding melt (taking apart finite chain-size effects). The test
chains interact with each other directly and through the density fluctuations of the sur-
rounding melt. The scaling of their effective interaction may be obtained from the potential
of mean force U(r, s) ≡ − ln(g(r, s)/g(∞, s)) where g(r) is the probability to find the second
chain at a distance r assuming the first segment at the origin (r = 0). Since the correlation
hole is shallow for large s, expansion leads to U(r, s) ≈ 1− g(r, s)/g(∞, s) ≈ c(r, s)/ρ with
c(r, s) being the density distribution of a test chain around its center of mass. This distri-
bution scales as c(r ≈ 0, s) ≈ s/R(s)d close to the center of mass (d being the dimension of
space) and decays rapidly at distances of order R(s) [5]. Hence, the interaction strength at
r/R(s) ≪ 1 is set by u(s) ≡ U(0, s) ≡ c(0, s)/ρ ≈ s/ρR(s)d ∼ s1−d/2 [12, 19]. Interestingly,
u(s) does not depend explicitly on the bulk compression modulus v. It is dimensionless and
independent of the definition of the monomer unit, i.e. it does not change if λ monomers
are regrouped to form an effective monomer (ρ → ρ/λ, s → s/λ) while keeping the segment
size R fixed.
Connectivity and swelling. To connect both test chains to form a chain of length 2s
the effective energy u(s) has to be paid and this repulsion will push the half-segments
apart. We consider next a segment of length s in the middle of a very long chain. All
interactions between the test segment and the rest of the chain are first switched off but we
keep all other interactions, especially within the segment and between the segment monomers
and monomers of surrounding chains. The typical size R(s) of the test segment remains
essentially unchanged from the size of an independent chain of same strand length. If we
now switch on the interactions between the segment and monomers on adjacent segments
of same length s, this corresponds to an effective interaction of order u(s) as before. (The
effect of switching on the interaction to all other monomers of the chain is inessential at
scaling level, since these other monomers are more distant.) Since this repels the respective
segments from each other, the corresponding subchain is swollen compared to a Gaussian
chain of non-interacting segments. It is this effect we want to characterize.
Perturbation approach in three dimensions. In the following we will exclusively consider
chain segments s which are much larger than the number of monomers g ≡ 1/vρ contained
in a blob [5], i.e. we will look on a scale where incompressibility matters. (The number g
is also sometimes called “dimensionless compressibility” [14].) Interestingly, when taken at
s = g the interaction strength takes the value
u(s = g) ≈
ρbdeg
(vρ)d/2−1
≈ Gz (6)
with Gz being the standard Ginzburg parameter used for the perturbation calculation
of strongly interacting polymers [6]. Hence, the segmental correlation hole potential
u(s) ≈ Gz(g/s)d/2−1 ≪ Gz for d > 2 and s ≫ g. Although for real polymer melts as
for computational systems large values of Gz ≈ 1 may sometimes be found, u(s) ∼ 1/
decreases rapidly with s in three dimensions, as illustrated in Fig. 2(b), and standard per-
turbation calculations can be successfully performed.
As sketched in the next paragraph these calculations yield quantities K[u] which are
defined such that they vanish (K[u = 0] = 0) if the perturbation potential u(s) is switched
off and are then shown to scale, to leading order, linearly with u. For instance, for the
quantity Kp(s), defined in Eq. (1), characterizing the deviation of the chain segment size
from Flory’s hypothesis one thus expects the scaling
Kp[u(s)] ≈ +u(s) ≈ +
ρRd(s)
. (7)
The +-sign indicated marks the fact that the prefactor has to be positive to be consis-
tent with the expected swelling of the chains. Consequently, the typical segment size,
R(s)/be
s ≈ 1 − u(s), must approach the asymptotic limit for large s from below. For
three dimensional solutions Eq. (7) implies that Kp(s) should vanish rapidly as 1/(ρb
(This is different in thin films where u(s) ≈ Gz decays only logarithmically [12] as may be
seen from Eq. (33) given below.) Taking apart the prefactors — which require a full calcu-
lation — this corresponds exactly to Eq. (5) with a swelling coefficient ce ≈ cp ≈ 1/ρb3e in
agreement with Eq. (4). Note also that the predicted deviations are inversely proportional
to b3e , i.e. the more flexible the chains, the more pronounced the effect. Similar relations
K[u] ∼ u may also be formulated for other quantities and will be tested numerically in
Sec. IV. There, we will also check that the linear order is sufficient.
B. Perturbation calculation
Generalities. Before delving more into our computational results we summarize here how
Eqs. (3-5) and related relations have been obtained using standard one-loop perturbation
calculation. The general task is to determine 〈A〉 ≈ 〈A〉0 (1+ 〈U〉0)−〈AU〉0 for measurable
quantities A such as the squared distance between two monomers n andm on the same chain,
A = r2nm. Here, 〈...〉0 denotes the average over the distribution function of the unperturbed
ideal chain of bond length b and U =
dl ṽ(rkl) the effective perturbation potential.
We discuss first the general results in the scale free regime (1 ≪ s ≪ N), argue then that b
should be renormalized to the effective bond length be and sketch finally the calculation of
finite chain-size effects.
The scale free regime. Following Edwards [6, 7, 8], the Gaussian (or “Random Phase”
[5]) approximation of the pair interaction potential in real space is
ṽ(r) = v
δ(r)− exp(−r/ξ)
4πrξ2
where v is a parameter which tunes the monomer interaction. (It is commonly associated
with the bare excluded volume of the monomers [6], but should more correctly be identified
with the bulk modulus effectively measured for the system. See the discussion of Eq. (15) of
Ref. ([13]).) The effective potential consists of a strongly repulsive part vδ(r) of very short
range, and a weak attractive part of range ξ where the correlation length of the density
fluctuations is given by ξ2 = b2g/12 with g = 1/ρv. In Fourier space Eq. (8) is equivalent to
ṽ(q) = v
q2 + ξ−2
with q being the wave vector. This is sufficient for calculating the scale free regime corre-
sponding to asymptotically long chains where chain end effects may be ignored. The different
graphs one has to compute are indicated in Fig. 1. For A = r2 (with 1 ≪ n < m ≪ N) this
yields, e.g.,
(m− n)−
24/π3
= b2es
. (10)
In the second line we have used the definition of the swelling coefficient ce indicated in
Eq. (4) and have set
b2e ≡ b2
1 + p
1 + p
with Gz ≡
vρ/b3ρ and p = 1. (The prefactor p has been added for convenience.) The
coefficient be of the leading Gaussian term in Eq. (10) — entirely due to the graph Ii
describing the interactions of monomers inside the segment between n and m — has been
predicted long ago by Edwards [6]. It describes how the effective bond length is increased
from b to be under the influence of a small excluded volume interaction. The second term in
Eq. (10) entails the 1/
s-swelling which is investigated numerically in this paper. It does
only depend on b and ρ but, more importantly, not on v — in agreement with the scaling
of u(s) discussed in Sec. IIA. The relative weights contributing to this term are indicated
in Fig. 1 in units of −
6/π3vξ2/b3
s. The diagrams I− and I+ are obviously identical in
the scale free limit. Note that the interactions described by the strongest graph Ii align the
bonds ln and lm while the others tend to reduce the effect.
For higher moments of the segment size distribution G(r, s) it is convenient to calculate
first the deviation of the Fourier-Laplace transformation of δG(r, s) and to obtain the mo-
ments from the coefficients of the expansion of this “generating function” in terms of the
squared wave vector q2. As explained in detail in the Appendix A this yields more generally
(2p+ 1)!
(b2es)
3(2pp!p)2
2(2p+ 1)!
)2p−3
where we have used Eq. (11) with general p. Obviously, Eq. (12) is consistent with our
previous finding Eq. (10) for p = 1. The corresponding segmental size distribution is
G(r, s) =
2πb2es
2πb2s
f(n) (13)
with n = r/b
s and f(n) being the same function as indicated in Eq. (3). The leading
Gaussian terms in Eqs. (12) and (13) depend on the effective bond length be, the second
only on the Kuhn length b of the reference chain. When comparing these result with Eqs. (3)
and (5) proposed in the Introduction, one sees that both equations are essentially identical —
taken apart, however, that they depend on b and be. Note the conspicuous factor (b/be)
in Eq. (12) which would strongly reduce the empirical swelling coefficients cp = ce(b/be)
for large p if b and be were different.
Interpretation of first-loop results in different contexts. The above perturbation results
may be used directly to describe the effect of a weak excluded volume v on a reference system
of perfectly ideal polymer melts with Kuhn segment length b where all interactions have been
switched off (v = 0). It is expected to give a good estimation for the effective bond length
be only for a small Ginzburg parameter: Gz ≪ 1. For the dense melts we want to describe
this does not hold (Sec. III) and one cannot hope to find a good quantitative agreement
with Eq. (11). Note also that large wave vectors contribute strongly to the leading Gaussian
term. The effective bond length be is, hence, strongly influenced by local and non-universal
effects and is very difficult to predict in general.
Our much more modest goal is to predict the coefficient of the 1/
s-perturbation and
to express it in terms of a suitable variational reference Hamiltonian characterized by a
conveniently chosen Kuhn segment b and the measured effective bond length be (instead of
Eq. (11)). Following Muthukumar and Edwards [10], we argue that for dense melts b should
be renormalized to be to take into account higher order graphs. No strict mathematical
proof can be given at present that the infinite number of possible graphs must add up in
this manner. Our hypothesis relies on three observations:
• The general scaling argument discussed in Sec. IIA states that we have only one
relevant length scale in this problem, the typical segment size R(s) ≈ be
s itself. The
incompressibility constraint cannot generate an additional scale. It is this size R(s)
which sets the strength of the effective interaction which then in turn feeds back to
the deviations of R(s) from Gaussianity. Having a bond length b in addition to the
effective bond length be associated with R(s) would imply incorrectly a second length
scale b
s varying independently with the bulk modulus v. (We will check explicitly
below in Fig. 13 that there is only one length scale.) This implies b/be = const v
• Thus, since by construction b/be = 1 for v → 0, it follows that both lengths should be
equal for all v.
• We know from Eq. (12) that the empirical coefficients cp = ce(b/be)2p−3 should depend
strongly on the moment considered if the ratio b/be is not close to unity. It will be
shown below (Fig. 6) that cp/ce ≈ 1 for all p. This implies b ≈ be.
Finite chain size effects. To describe properly finite chain size corrections Eq. (9) must
be replaced by the general linear response formula
ṽ(q)
+ ρF (q) (14)
with F (q) = NfD(x) being the form factor of the Gaussian reference chain given by Debye’s
function fD(x) = 2(e
−x − 1 + x)/x2 with x = (qb)2N/6 [6]. This approximation allows
in principle to compute, for instance, the (mean-squared) total chain end-to-end distance,
A = (rN −r1)2. One verifies readily (see [6], Eq. (5.III.9)) that the effect of the perturbation
may be expressed as
〈A〉0 〈U〉0 − 〈AU〉0 =
(2π)3
ṽ(q)
ds s2(N − s) exp
q2b2s
. (15)
We take now first the integral over s. In the remaining integral over q small q wave vectors
contribute to the
N -swelling while large q renormalize the effective bond length of the
dominant Gaussian behaviour linear in N (as discussed above). Since we wish to determine
the non-Gaussian corrections, we may focus on small wave vectors q ≪ 1/ξ. Since in this
limit 1/v = ρg ≪ ρF (q), one can neglect in Eq. (14) the 1/v contribution to the inverse
effective interaction potential. We thus continue the calculation using the much simpler
ṽ(N, x) = 1/(NρfD(x)). This allows us to express the swelling as
〈(rN − r1)2〉
I(xu). (16)
To simplify the notation we have set here finally b = be in agreement with the hypothesis
discussed above. The numerical integral I(xu) =
dx . . . over x is slowly convergent at
infinity. As a consequence the estimate I(∞) = 1.59 may be too large for moderate chain
lengths. In practice, convergence is not achieved for values xu(N) ≈ (b/ξ)2N corresponding
to the screening length ξ.
We remark finally that numerical integration can be avoided for various properties if the
Padé approximation of the form factor, F (q) = N/(1 + (qb)2N/12), is used. This allows
analytical calculations by means of the simplified effective interaction potential
ṽ(q) =
12ρb3
. (17)
This has been used for instance for the calculation of finite chain size effects for the bond-
bond correlation function discussed in Sec. IVC below [61].
III. COMPUTATIONAL MODELS AND TECHNICAL DETAILS
A. Bond fluctuation model
A widely-used lattice Monte Carlo scheme for coarse-grained polymers. The body of
our numerical data comes from the three dimensional bond fluctuation model (BFM) —
a lattice Monte Carlo (MC) algorithm where each monomer occupies eight sites of a unit
cell of a simple cubic lattice [24, 25, 26]. Our version of the BFM with 108 bond vectors
corresponds to flexible athermal chain configurations [23]. All length scales are given in units
of the lattice constant and time in units of Monte Carlo Steps (MCS). We use cubic periodic
simulation boxes of linear size L = 256 containing nmon = ρL
3 = 220 ≈ 106 monomers.
This monomer number corresponds to a monomer number density ρ = 0.5/8 where half of
the lattice sites are occupied (volume fraction 0.5). The large system sizes used allow us to
suppress finite box-size effects for systems with large chains. Using a mix of local, slithering
snake [27, 28, 29], and double-bridging [17, 23, 30, 31] MC moves we were able to equilibrate
dense systems with chain lengths up to N = 8192.
Equilibration and sampling of high-molecular BFM melts. Standard BFM implementa-
tions [26, 32, 33] use local MC jumps to the 6 closest lattice sites to prevent the crossing
of chains and conserve therefore the chain topology. These “L06” moves lead to very large
relaxation times, scaling at least as τe ∼ N3, as may be seen from Fig. 3 (stars). The relax-
ation time τe = R
e/6Ds indicated in this figure has been estimated from the self-diffusion
coefficient Ds obtained from the mean-square displacements of all monomers in the free
diffusion limit. (For the largest chain indicated for L06 dynamics only a lower bound for
τe is given.) Instead of this more realistic but very slow dynamical scheme we make jump
attempts to the 26 sites of the cube surrounding the current monomer position (called “L26”
moves). This allows the chains to cross each other which dramatically speeds up the dynam-
ics, especially for long chains (N > 512). If only local moves are considered, the dynamics
is perfectly consistent with the Rouse model [6]. As shown in Fig. 3, we find τe ≈ 530N2 for
L26 dynamics. This is, however, still prohibitive by large for sampling configurations with
the longest chain length N we aim to characterize [62].
Slithering snake moves. In addition to the local moves one slithering snake move per
chain is attempted on average per MCS corresponding to the displacement of N monomer
along the chain backbone. Note that in our units two spatial displacement attempts per
MCS are performed on average per monomer, one for a local move and one for a snake
move. (In practice, it is computationally more efficient for large N to take off a monomer
at one chain end and to paste it at the other leaving all other monomers unaltered. Before
dynamical measurements are performed the original order of beads must then be restored.)
Interestingly, a significantly larger slithering snake attempt frequency would not be useful
since the relaxation time of slithering snakes without or only few local moves increases
exponentially with mass [29, 34] due to the correlated motion of snakes [35]. In order to
obtain an efficient free snake diffusion (with a chain length independent curvilinear diffusion
coefficient Dc(N) ∼ N0 and τe ≈ N2/Dc(N) ∼ N2 [28, 29]) it is important to relax density
fluctuations rapidly by local dynamical pathways. As shown in Fig. 3 (squares), we find
a much reduced relaxation time τe ≈ 40N2 which is, however, still unconveniently large
for our longest chains. Note that most of the CPU time is still used by local moves. The
computational load per MCS remains therefore essentially chain length independent.
Advantages and pitfalls of double-bridging moves. Double-bridging (DB) moves are very
useful for high densities and help us to extend the accessible molecular masses close to
104. As for slithering snake moves we use all 108 bond vectors to switch chain segments
between two different chains. Only chain segments of equal length are swapped to conserve
monodispersity. Topolocial constraints are again systematically and deliberately violated.
Since more than one swap partner is possible for a selected first monomer, delicate detailed
balance questions arise. This is particularly important for short chains and is discussed
in detail in Ref. [23]. Technically, the simplest solution to this problem is to refuse all
moves with more than one swap partner (to be checked both for forward and back move).
The configurations are screened with a frequency fDB for possible DB moves where we
scan in random order over the monomers. The frequency should not be too large to avoid
(more or less) immediate back swaps and monomers should move at least out of the local
monomer cage and over a couple of lattice sites. We use fDB = 0.1 between DB updates for
the configurations reported here. (The influence of fDB on the performance has not been
explored systematically, but preliminary results suggest a slightly smaller DB frequency for
future studies.) The diffusion times over the end-to-end distance for this case are indicated
in Tab. I. As shown in Fig. 3, we find empirically τe(N) ≈ 13N1.62. For N = 8192 this
corresponds to 3 · 107 MCS. This allows us even for the largest chain lengths to observe
monomer diffusion over several Re within the 10
8 MCS which are feasible on our XEON-PC
processor cluster.
The efficiency of DB moves is commonly characterized in terms of the relaxation time
τee of the end-to-end vector correlation function [30, 31]. For normal chain dynamics this
would indeed characterize the longest relaxation time of the system, i.e. τe ≈ τee. For the
double-bridging this is, however, not sufficient since density fluctuations do not couple to the
bridging moves and can not be relaxed. We find therefore that configurations equilibrate
on time scales given by τe rather than by τee ≪ τe. This may be verified, for instance,
from the time needed for the distribution Re(s) (and especially its spatial components) to
equilibrate. The criterion given in the literature [30] is clearly not satisfactory and may lead
to insufficiently equilibrated configurations. In summary, equilibration with DB moves still
requires monomer diffusion over the typical chain size, however at a much reduced price.
Some properties of our configurations. The Tables I and II summarize some system
properties obtained for our reference density ρ = 0.5/8. Averages are performed over all
chains and 1000 configurations. These configurations may be considered to be independent
for N < 4096. Only a few independent configurations exist for the largest chain length
N = 8192 which has to be considered with some care. Taking apart this system, chains are
always much smaller than the box size. For asymptotically long chains, we obtain an average
bond length 〈|l|〉 ≈ 2.604, a root-mean-squared bond length l ≡ 〈l2〉1/2 ≈ 2.635 and an
effective bond length be ≈ 3.244 — as we will determine below in Sec. IVA. This corresponds
to a ratio C∞ ≡ b2e/l2 ≈ 1.52 and, hence, to a persistence length lp = l(C∞ + 1)/2 ≈ 3.32
[23]. Especially, we find from the zero wave vector limit of the total structure factor S(q)
a low (dimensionless) compressibility g = S(q → 0)/ρ ≈ 0.246 which compares well with
real experimental melts. From the measured bulk compression modulus v ≡ 1/g(ρ)ρ and
the effective bond length be one may estimate a Ginzburg parameter Gz =
vρ/b3eρ ≈ 0.96.
Following Ref. [13] the interaction parameter v is supposed here to be given by the full
inverse compressibility and not just by the second virial coefficient.
B. Bead spring model
Hamiltonian. Additionally, molecular dynamics simulations of a bead-spring model
(BSM) [36] were performed to dispel concerns that our results are influenced by the un-
derlying lattice structure of the BFM. The model is derived from a coarse-grained model
for polyvinylalcohol which has been employed to study polymer crystallization [37]. It is
characterized by two potentials: a non-bonded potential of Lennard-Jones (LJ) type and a
harmonic bond potential. While the often employed Kremer-Grest model [38] uses a 12− 6
LJ potential to describe the non-bonded interactions Unb(r), our non-bonded potential has
a softer repulsive part. It is given by
Unb(r) = 1.511
, (18)
which is truncated and shifted at the minimum at rmin ≈ 1.15. Note that all length scales
are given in units of σ0 and we use LJ units [39] for all BSM data (mass m = 1, Boltzmann
constant kB = 1). The parameters of the bond potential, Ub(r) = 1120(r− lb)2, are adjusted
so that the average bond length l(ρ = 0.84) ≈ lb = 0.97 is approximately the same as in the
standard Kremer-Grest model [38]. The average bond length and the root-mean-squared
bond length are almost identical for the BSM due to the very stiff bond potential. Since
rmin/l ≈ 1.16 bonded monomers penetrate each other significantly.
Equilibration and sampling. We perform standard molecular dynamics simulations in the
canonical ensemble with a Langevin thermostat (friction constant Γ = 0.5) at temperature
T = 1. The equations of motion are integrated by the velocity-Verlet algorithm [39]. To
improve the statistics for large chain length, we have implemented additional double-bridging
moves. Since only few of these MC moves are accepted per unit time, this does affect neither
the stability nor the accuracy of the molecular dynamics sweeps.
Some properties obtained. For clarity, we show only data for chain length N = 1024
and number density ρ = 0.84, the typical melt density of the Kremer-Grest model [38].
For the reported data we use periodic simulation boxes of linear size L ≈ 62 containing
nmon = 196608 monomers, but we have also sampled different boxes sizes (up to L = 77.5)
to check for finite box-size effects. For the reference density a dimensionless compressibility
g ≈ 0.08 is found which is about three times smaller than for our BFM melt. For the
effective bond length we obtain be ≈ 1.34, i.e. BSM chains (C∞ ≈ 1.91, lp ≈ 1.41) are
slightly stiffer than the corresponding BFM polymers. Fortunately, the product ρb3e ≈ 2 is
roughly similar in both models and one expects from Eq. (4) a similar swelling for large
s. Note finally that the Ginzburg parameter Gz ≈ 1.8 is much larger than for the BFM
systems. As we have emphasized in Sec. II, this should, however, not influence the validity
of the perturbation prediction of the expected 1/
s-swelling of the chains when expressed
in terms of the measured effective bond length.
IV. NUMERICAL RESULTS
As illustrated in Fig. 1, a chain segment of curvilinear length s > 0 is identified by two
monomers n and m = n + s on the same chain. We compute here various moments of
chain segment properties where we ensemble-average over all chains and all start points n.
The statistical accuracy must therefore always decrease for large s. We concentrate first on
the second moment (p = 1) of the segmental size distribution. Higher moments and the
segmental size distribution are discussed in Sec. IVE.
A. The swelling of chain segments
Scale free regime for 1 ≪ s ≪ N . The mean-squared segment size R2e(s) = 〈r2〉 is
presented in the Figs. 4, 5 and 6. The first plot shows clearly that chain segments are swollen,
i.e. R2e(s)/s increases systematically and this up to very large curvilinear distances s. Only
BFM data are shown for clarity. A similar plots exists for the BSM data. In agreement
with Eq. (5) for p = 1, the asymptotic Gaussian behavior (dashed line) is approached
from below and the deviation decays as u(s) ∝ 1/
s (bold line). The bold line indicated
corresponds to be = 3.244 and c1 ≈ ce ≈ 0.41 which fits nicely the data over several decades
in s — provided that chain end effects can be neglected (s ≪ N). Note that a systematic
underestimation of the true effective bond length would be obtained by taking simply the
largest R2e(s)/s ≈ 3.232 value available, say, for monodisperse chains of length N = 2048.
Finite chain-size effects. Interestingly, R2e(s)/s does not approach the asymptotic limit
monotonicly. Especially for short chains one finds a non-monotonic behavior for s → N .
This means that the total chain end-to-end distance Re(s = N − 1) must show even more
pronounced deviations from the asymptotic limit. This is confirmed by the dashed line
representing the b2e(N) ≡ R2e(N − 1)/(N − 1) data points given in Tab. I. We emphasize
that the non-monotonicity of R2e(s)/s becomes weaker with increasing N and that, as one
expects, the inner distances, as well as the total chain size, are characterized by the same
effective bond length be for large s or N . The non-monotonic behavior may be qualitatively
understood by the reduced self-interactions at chain ends which lessens the swelling on these
scales. These finite-N corrections have been calculated analytically using the full Debye
function for the effective interaction potential ṽ(q), Eq. (14). The prediction for the total
chain end-to-end vector given in Eq. (16) is indicated in Fig. 4 (dash-dotted line) where we
have replaced the weakly N -dependent integral I(xu) by its upper bound value for infinite
chains
R2e(N − 1)
b2e (N − 1)
1.59ce√
N − 1
. (19)
We have changed here the chain length N in the analytical formula (obtained for large chains
where N ≈ N − 1) to the curvilinear length N − 1. This is physically reasonable and allows
to take better into account the behavior of small chains. Note that Eq. (19) is similar to
Eq. (5) — apart from a slightly larger prefactor explaining the observed stronger deviations.
Theory compares well with the measured data for large N . It does less so for smaller N ,
as expected, where the chain length dependence of the numerical integral I(xu(N)) ≤ 1.59
must become visible. This explains why the data points are above the dash-dotted line.
Note also that additional non-universal finite-N effects not accounted for by the theory are
likely for small N . In contrast to this, Re(s) is well described by the theory even for rather
small s provided that N is large and chain end effects can be neglected. In summary, it is
clear that one should use the segment size Re(s) rather than the total chain size to obtain
in a computational study a reliable fit of the effective bond length be.
Extrapolation of the effective bond length of asymptotically long chains. The represen-
tation chosen in Fig. 4 is not the most convenient one for an accurate determination of be
and c1. How precise coefficients may be obtained according to Eq. (5) is addressed in the
Figs. 5 and 6. The fitting of the effective bond length be and its accuracy is illustrated in
Fig. 5 for BFM chains of length N = 2048. This may be first done approximately in linear
coordinates by plotting R2e(s)/s as a function of 1/
s (not shown). Since data for large s
are less visible in this representation, we recommend for the fine-tuning of be to switch then
to logarithmic coordinates with a vertical axis y = 1 − R2e(s)/b2es for different trial values
of be. The correct value of be is found by adjusting the vertical axis y such that the data
extrapolates linearly as a function of 1/
s to zero for large s. We assume for the fine-tuning
that higher order perturbation corrections may be neglected, i.e. we take Eq. (5) literally.
(We show below that higher order corrections must indeed be very small.) The plot shows
that this method is very sensitive, yielding a best value that agrees with the theory over
more than one order of magnitude without curvature. As expected, it is not possible to
rationalize the numerically obtained values be ≈ 3.244 for the BFM and be ≈ 1.34 for the
BSM using Eq. (11). According to Eq. (4) these fit values imply the theoretical swelling
coefficients ce = 0.41 for the BFM and ce = 0.44 for the BSM.
Empirical swelling coefficients. As a next step the horizontal axis is rescaled such that
all data sets collapse on the bisection line, i.e. using Eq. (5) we fit for the empirical swelling
coefficient c1 and compare it to the predicted value ce. This rescaling of the axes allows to
compare both models in Fig. 6. For clarity the BSM data have been shifted upwards. For
the BFM we find c1/ce ≈ 1.0, as expected, while our BSM simulations yield a slightly more
pronounced swelling with c1/ce ≈ 1.2.
Segmental radius of gyration. Also indicated in Fig. 6 is the segmental radius of gyration
Rg(s) (filled circles) computed as usual [6] as the variance of the positions of the segment
monomers around their center of mass. Being the sum over all s + 1 monomers, it has a
much better statistics compared to Re(s). The scaling used can be understood by expressing
the radius of gyration R2g(s) =
(s+1)2
l=n r
in terms of displacement vectors rkl
[6]. Using Eq. (5) and integrating twice this yields
6R2g(s)
b2e(s+ 1)
. (20)
Plotting the l.h.s. of this relation against the r.h.s. we obtain a perfect data collapse on the
bisection line where we have used the same parameters be and c1 as for the mean-squared
segment size. This is an important cross-check which we strongly recommend. Different
values indicate insufficient sample equilibration.
B. Chain connectivity and recursion relation
As was emphasized in Sec. IIA the observed swelling is due to an entropic repulsion
between chain segments induced by the incompressibility of the melt. To stress the role of
chain connectivity we repeat the general scaling argument given above in a form originally
proposed by Semenov and Johner for ultrathin films [12]. As shown in Fig. 7 we test the
relation
Kλ(s) ≡
R2e(λs)− λR2e(s)
λ)R2e(s)
≈ u(s) ≡ s
ρRe(s)d
with Kλ(s) being a direct measure of the non-Gaussianity (λ being a positive number)
comparing the size of a segment of length λs with the size of λ segments of length s joined
together. (The prefactor 1/(λ −
λ) in the definition of Kλ(s) has been introduced for
convenience.) Equivalently, this can be read as a measure for the swelling of a chain where
initially the interaction energy u between the segments has been switched off. Kλ is a
functional of u(s) with Kλ[u = 0] = 0. The analytic expansion of the functional must be
dominated by the linear term (as indicated by ≈ in the above relation) simply because u is
very small. Altogether, Eq. (21) yields a recursion relation relating Re(λs) with Re(s) for
any λ provided 1 ≪ s < λs ≪ N . It can be solved, leading (in lowest order) to Eq. (5) with
p = 1. This may be seen from the ansatz R2e(s) = b
es(1− ce/sω−1+ . . .) which readily yields
ω = 3/2 and ce ≈ 1/ρb3e .
Eq. (21) has been validated directly in Fig. 7 for λ = 2 (corresponding to two segments
of length s joined together) for the BFM and the BSM as indicated. In addition, for BFM
chains of length N = 2048 several values of λ have been given. As suggested by Eq. (4),
we have plotted Kλ(s) as a function of (c1/ce) u(s) with u(s) ≡
24/π3s/ρR3e(s) ≈ ce/
The prefactor of u(s) allows a convenient comparison with Fig. 5. Note the perfect data
collapse for all data sets. More importantly, the predicted linearity is well confirmed for
large segments (1 ≪ s) and this without any tunable parameter for the vertical axis, as was
needed in the previous Figs. 5 and 6.
C. Intrachain bond-bond correlations
Expectation from Flory’s hypothesis. An even more striking violation of Flory’s ideality
hypothesis may be obtained by computing the bond-bond correlation function, defined by
the first Legendre polynomial P (s) = 〈lm=n+s · ln〉 /l2 where the average is performed, as
before, over all possible pairs of monomers (n,m = n+ s) [63]. Here, li = ri+1 − ri denotes
the bond vector between two adjacent monomers i and i + 1 and l2 = 〈l2n〉n the mean-
squared bond length. The bond-bond correlation function is generally believed to decrease
exponentially [4]. This belief is based on the few simple single chain models which have
been solved rigorously [4, 40] and on the assumption that all long range interactions are
negligible on distances larger than the screening length ξ. Hence, only correlations along
the backbone of the chains are expected to matter and it is then straightforward to work
out that an exponential cut-off is inevitable due to the multiplicative loss of any information
transferred recursively along the chain [4].
Asymptotic behavior in the melt. That this reasoning must be incorrect follows imme-
diately from the relation
P (s) =
R2e(s) (22)
expressing the bond-bond correlation function as the curvature of the second moment of
the segment size distribution. It is obtained from the identity 〈ln · lm〉 = 〈∂nrn · ∂mrm〉 =
−∂n∂m 〈r2nm〉 /2. (Note that the velocity correlation function is similarly related to the
second derivative of the mean-square displacement with respect to time [41].) Hence, P (s)
allows us to probe directly the non-Gaussian corrections without any ideal contribution.
This relation together with Eq. (5) suggests an algebraical decay P (s) = cP/s
ω with
ω = 3/2 , cP = c1 (be/l)
2/8 ≈
ρl2be
of the bond-bond correlation function for dense solutions and melts, rather than the ex-
ponential cut-off expected from Flory’s hypothesis. This prediction (bold line) is perfectly
confirmed by the larger chains (N > 256) indicated in Fig. 8. In principle, the swelling
coefficient, c1 ∼ cP, may also be obtained from the power law amplitude of the bond-bond
correlation function, however, to lesser accuracy than by the previous method (Fig. 5). One
reason is that P (s) decays very rapidly and does not allow a precise fit beyond s ≈ 102. The
values of cP obtained from c1 are indicated in Tab. II. Data from the BSM have also been
included in the figure to demonstrate the universality of the result. The vertical axis has
been rescaled with cP which allows to collapse the data of both models.
Finite chain-size corrections. As can be seen for N = 16, exponentials are compatible
with the data of short chains. This might explain how the power law scaling has been
overlooked in previous numerical studies, since good statistics for large chains (N > 1000)
has only become available recently. However, it is clearly shown that P (s) approaches
systematically the scale free asymptote with increasing N . The departure from this limit is
fully accounted for by the theory if chain end effects are carefully considered (dashed lines).
Generalizing Eq. (23) and using the Padé approximation, Eq. (17), perturbation theory
yields
P (s) =
1 + 3x+ 5x2
1 + x
(1− x)2 (24)
where we have set x =
s/N . For x ≪ 1 this is consistent with Eq. (23). In the limit of
large s → N , the correlation functions vanish rigorously as P (s) ∝ (1 − x)2. Considering
that non-universal features cannot be neglected for short chain properties and that the
theory does not allow for any free fitting parameter, the agreement found in Fig. 8 is rather
satisfactory.
D. Higher moments and associated coefficients
Effective bond length and empirical swelling coefficients. The preceding discussion fo-
cused on the second moment of the segmental size distribution G(r, s). We have also com-
puted for both models higher moments 〈r2p〉 with p ≤ 5. If traced in log-linear coordinates
as y = (6pp! 〈r2p〉 /(2p + 1)!sp)1/p vs. x = s higher moments approach b2e from below —
just as the second moment presented in Fig. 4. The deviations from ideality are now more
pronounced and increase with p (not shown). The moments are compared in Fig. 6 with
Eq. (5) where they are rescaled as y = Kp(s) as defined in Eq. (1) and plotted as functions
of x =
3(2pp!p)2
2(2p+1)!
. The prediction is indicated by bold lines. It is important that the same
effective bond length be is obtained from the analysis of all functions Kp(s) as illustrated in
Fig. 6. Otherwise we would regard equilibration and statistics as insufficient.
The empirical swelling coefficients cp are obtained, as above in Sec. IVA, by shifting the
data horizontally. A good agreement with the expected cp/ce ≈ 1 is found for both models
and all moments as may be seen from Tab. II. This confirms the renormalization of the
Kuhn segment b → be of the Gaussian reference chain in agreement with our discussion in
Sec. II B. Otherwise we would have measured empirical coefficients decreasing strongly as
cp/ce ≈ (b/be)2p−3 with p. Since the effective bond length of non-interacting chains are known
for the BFM (b ≈ 2.688) and the BSM (b ≈ 0.97), one can simply check, say for p = 5, that
the non-renormalized values would correspond to the ratios c5/ce ≈ (2.688/3.244)7 ≈ 0.3 for
the BFM and c5/ce ≈ (0.97/1.34)7 ≈ 0.1 for the BSM. This is clearly not consistent with
our data.
It should be emphasized that both coefficients be and cp are more difficult to determine
for large p, since the linear regime for x ≪ 1 in the representation chosen in Fig. 6 becomes
reduced. For large x ≫ 1 one finds that y(x) → 1, i.e. 〈r2p〉 /b2pe sp → 0. This trivial
departure from both Gaussianity and the 1/
s-deviations we try to describe, is due to the
finite extensibility of chain segments of length s which becomes more marked for larger
moments probing larger segment sizes. The data collapse for both x-regimes is remarkable,
however. Incidentally, it should be noted that for the BSM the empirical swelling coefficients
are slightly larger than expected. At present we do not have a satisfactory explanation for
this altogether minor effect, but it might be attributed to the fact that neighbouring BSM
beads along the chain strongly interpenetrate — an effect not considered by the theory.
Non-Gaussian parameter αp. The failure of Flory’s hypothesis can also be demonstrated
by means of the standard non-Gaussian parameter
αp(s) ≡ 1−
(2p+ 1)!
〈r2pnm〉
〈r2nm〉
p (25)
comparing the 2p-th moment with the second moment (p = 1). In contrast to the closely
related parameter Kp(s) this has the advantage that here two measured properties are com-
pared without any tuneable parameter, such as be, which has to be fitted first. Fig. 9 presents
αp(s) vs. ce/
s for the three moments with p = 2, 3, 4. For each p we find perfect data col-
lapse for all chain lengths and both models and confirm the linear relationship αp(s) ≈ u(s)
expected. The lines indicate the theoretical prediction
αp(s) =
3 (2pp!p)2
2 (2p+ 1)!
which can be derived from Eq. (5) by expanding the second moment in the denominator. An
alternative derivation based on the coefficients of the expansion of the generating function
G(q, s) in q2 is indicated by Eq. (A2) in the Appendix. Having confirmed above that cp/ce ≈
1, we assume in Eq. (26) that cp = ce to simplify the notation. The prefactors 6/5, 111/35
and 604/105 for p = 2, 3 and 4 respectively are nicely confirmed. They increase strongly
with p, i.e. the non-Gaussianity becomes more pronounced for larger moments as already
mentioned. Note also the curvature of the data at small s due to the finite extensibility
of the segments which becomes more marked for higher moments. If one plots αp(s) as a
function of the r.h.s. of Eq. (26) all data points for all moments and even for too small s
collapse on one master curve (not shown) — just as we have seen before in Fig. (6).
Correlations of different directions. A similar correlation function is presented in Fig. 10
which measures the non-Gaussian correlations of different spatial directions. It is defined by
Kxy(s) ≡ 1−
〈x2 y2〉
〈x2〉 〈y2〉
for the two spatial components x and y of the vector r as illustrated by the sketch given
at the bottom of Fig. 10. Symmetry allows to average over the three pairs of directions
(x, y), (x, z) and (x, z). Following the general scaling argument given in Sec. II we expect
Kxy(s) ≈ u(s) ≈ ce/
s which is confirmed by the perturbation result
Kxy(s) =
= K2(s). (28)
This is nicely confirmed by the linear relationship found (bold line) on which all data from
both simulation models collapse perfectly. The different directions of chain segments are
therefore coupled. As explained in the Appendix (Eq. (A3)), Kxy(s) and α2(s) must be
identical if the Fourier transformed segmental size distribution G(q, s) can be expanded
in terms of q2 and this irrespective of the values the expansion coefficients take. Fig. 10
confirms, hence, that our computational systems are perfectly isotropic and tests the validity
of the general analytical expansion.
The correlation function Kxy is of particular interest since the zero-shear viscosity should
be proportional to
∼ 〈x2y2〉 = 〈x2〉 〈y2〉 (1 − Kxy(s)). We assume here following
Edwards [6] that only intrachain stresses contribute to the shear stress σxy. Hence, our
results suggest that the classical calculations [6] — assuming incorrectly Kxy = 0 — should
be revisited.
E. The segmental size distribution
We turn finally to the segmental size distribution G(r, s) itself which is presented in
Figs. 11, 12 and 13. From the theoretical point of view G(r, s) is the most fundamental
property from which all others can be derived. It is presented last since it is computationally
more demanding — at least if high accuracy is needed — and coefficients such as be may
be best determined directly from the moments. The normalized histograms G(r, s) are
computed by counting the number of segment vectors between r − dr/2 and r + dr/2 with
dr being the width of the bin and one divides then by the spherical bin volume. Since
the BFM model is a lattice model, this volume is not 4πr2dr but given by the number of
lattice sites the segment vector can actually point to for being allocated to the bin. Incorrect
histograms are obtained for small r if this is not taken into account. (Averages are taken
over all segments and chains, just as before.) Clearly, non-universal physics must show up
for small vector length r and small curvilinear distance s and we concentrate therefore on
values r ≫ σ and s ≥ 31.
When plotted in linear coordinates as in Fig. 11, G(r, s) compares roughly with the
Gaussian prediction G0(r, s) given by Eq. (2), but presents a distinct depletion for small
segment sizes with n ≡ r/be
s ≪ 1 and an enhanced regime for n ≈ 1. A second depletion
region for large n ≫ 1 — expected from the finite extensibility of the segments — can be best
seen in the log-log representation of the data (not shown). To analyse the data it is better
to consider instead of G(r, s) the relative deviation δG(r, s)/G0(r, s) = G(r, s)/G0(r, s)− 1
which should further be divided by the strength of the segmental correlation hole, ce/
As presented in Fig. 12 this yields a direct test of the key relation Eq. (3) announced in the
Introduction. The figure demonstrates nicely the scaling of the data for all s and for both
models. It shows further a good collapse of the data close to the universal function f(n)
predicted by theory (bold line). Note that the depletion scales as 1/n for small segment
sizes (dashed line). The agreement of simulation and theory is by all standards remarkable.
(Obviously, error bars increase strongly for n ≫ 1 where G0(r, s) decreases strongly. The
regime for very large n where the finite extensibility of segments matters has been omitted
for clarity.) We emphasize that this scaling plot depends very strongly on the value be which
is used to calculate the Gaussian reference distribution.
If a precise value is not available we recommend to use instead the scaling variable
m = r/Re(s) for the horizontal axis, i.e. to replace the scale be
s estimated from the
behaviour of asymptotically long chains by the measured (mean-squared) segment size
for the given s. The Gaussian reference distribution is then accordingly G0(m,Re(s)) =
(3/2πRe(s)
2) exp(−3
m2). The corresponding scaling plot is given in Fig. 13. It is simi-
lar and of comparable quality as the previous plot. Changing the scaling variable from
n = r/be
s to m = r/Re(s) ≈ (r/be
s)(1 + ce/2
s) changes somewhat the universal func-
tion. Expanding the previous result, Eq. (3), this adds even powers of m to the function
f(n) given in Eq. (3)
f(n) ⇒ f(m) =
+ 9m+
m2 − 9
. (29)
That the two additional terms in the function are correct can be seen by computing the
second moment 4π
drr4δG(r, s) which must vanish by construction. The rescaled relative
deviation is somewhat broader than in the previous plot due to the additional term scaling
as m2. As already stressed this scaling does not rely on the effective bond length be and
is therefore more robust. It has the nice feature that it underlines that there is only one
characteristic length scale relevant for the swelling induced by the segmental correlation
hole, the typical size of the chain segment itself.
V. CONCLUSION
Issues covered and central theoretical claims. We have revisited Flory’s famous ideality
hypothesis for long polymers in the melt by analyzing both analytically and numerically
the segmental size distribution G(r, s) and its moments for chain segments of curvilinear
length s. We have first identified the general mechanism that gives rise to deviations from
ideal chain behavior in dense polymer solutions and melts (Sec. II). This mechanism rests
upon the interplay of chain connectivity and the incompressibility of the system which
generates an effective repulsion between chain segments (Fig. 2). This repulsion scales like
u(s) ≈ ce/
s where the “swelling coefficient” ce ≈ 1/b3eρ sets the strength of the interaction.
It is strong for small segment length s, but becomes weak for s → N in the large-N limit.
The overall size of a long chain thus remains almost ‘ideal’, whereas subchains are swollen as
described by Eq. (5). Most notably, the relative deviation δG(r, s)/G0(r, s) of the segmental
size distribution from Gaussianity should be proportional to u(s). As a function of segment
size r, the repulsion manifests itself by a strong 1/r-depletion at short distances r ≪ be
and a subsequent shift of the histogram to larger distances (Eq. (3)).
Summary of computational results. Using Monte Carlo and molecular dynamics simu-
lation of two coarse-grained polymer models we have verified numerically the theoretical
predictions for long and flexible polymers in the bulk. We have explicitly checked (e.g.,
Figs. 7, 9, 13) that the relative deviations from Flory’s hypothesis scale indeed as 1/
Especially, the measurement of the bond-bond correlation function P (s), being the second
derivative of the second moment of G(r, s) with respect of s, allows a very precise verifi-
cation (Fig. 8) and shows that higher order corrections beyond the first-order perturbation
approximation must be small. The most central and highly non-trivial numerical verification
concerns the data collapse presented in Figs. 12 and 13 for the segmental size distribution
of both computational models. All other statements made in this paper can be derived and
understood from this key finding. It shows especially that the swelling coefficient ce must
be close to the predicted value, Eq. (4).
It is well known [10] that the effective bond length is difficult to predict at low com-
pressibility and no attempt has been done to do so in this paper. We show instead how the
systematic swelling of chain segments – once understood – may be used to extrapolate for
the effective bond length of asymptotically long chains. Figs. 5 and 6 indicate how this may
be done using Eq. (5). The high precision of our data is demonstrated in Fig. 12 by the
successful scaling of the segmental size distribution.
For several moments 〈r2p〉 we have also fitted empirical swelling coefficients cp using
Eq. (5). In contrast to the effective bond length be these coefficients are rather well pre-
dicted by one-loop perturbation theory if the bond length b of the reference Hamiltonian is
renormalized to the effective bond length be, as we have conjectured in Sec. II B. Since the
empirical swelling coefficients, cp ≈ ce(b/be)2p−3, would otherwise strongly depend on the
moment taken, as shown in Eq. (12), our numerical data (Tab. II) clearly imply b/be ≈ 1.
Minor deviations found for the BSM samples may be attributed to the fact that monomers
along the BSM chains do strongly overlap — an effect not taken into account by the theory.
To clarify ultimately this issue we are currently performing a numerical study where we
systematically vary both the compressibility and the bond length of the BSM.
General background and outlook. The most striking result presented in this work con-
cerns the power law decay found for the bond-bond correlation function, P (s) ∝ 1/s3/2
(Fig. 8). This result suggests an analogy with the well-known long-range velocity correla-
tions found in dense fluids by Alder and Wainwright nearly fourty years ago [41, 42]. In
both cases, the ideal uncorrelated object is a random walker which is weakly perturbed (for
d > 2) by the self-interactions generated by global constraints. Although these constraints
are different (momentum conservation for the fluid, incompressibility for polymer melts) the
weight with which these constraints increase the stiffness of the random walker is always
proportional to the return probability. It can be shown that the correspondence of both
problems is mathematically rigorous if the fluid dynamics is described on the level of the
linearized Navier-Stokes equations [43].
We point out that the physical mechanism which has been sketched above is rather
general and should not be altered by details such as a finite persistence length — at least
not as long as nematic ordering remains negligible and the polymer chains are sufficiently
long. (Similarly, velocity correlations in dense liquids must show an analytical decay for
sufficiently large times irrespective of the particle mass and the local static structure of
the solution.) While this paper focused exclusively on scales beyond the correlation length
of the density fluctuations, i.e. qξ ≪ 1 or s/g ≫ 1, where the polymer solution appears
incompressible, effects of finite density and compressibility can be readily described within
the same theoretical framework and will be presented elsewhere [43]. To test our predictions,
flexible chains should be studied preferentially, since the chain length required for a clear-
cut description increases strongly with persistence length. This is in fact confirmed by
preliminary and on-going simulations using the BSM algorithm.
In this work we have only discussed properties in real space as a function of the curvilin-
ear distance s. These quantities are straightforward to compute in a computer simulation
but are barely experimentally relevant. The non-Gaussian deviations induced by the seg-
mental correlation hole arise, however, also for an experimentally accessible property, the
intramolecular form factor (single chain scattering function) F (q). As explained at the
end of the Appendix, the form factor can be readily obtained by integrating the Fourier
transformed segmental size distribution given in Eq. (3). This yields
q2F (q) ≈ 12
in agreement with the result obtained in Refs. [14, 19] by direct calculation of the form factor
for very long equilibrium polymers. As a consequence of this, the Kratky plot (q2F (q) vs.
wave vector q) should not exhibit the plateau expected for Gaussian chains in the scale-free
regime, but rather noticeable non-monotonic deviations. See Fig. 3 of [19]. This result
suggests to revisit experimentally this old pivotal problem of polymer science.
Our work is part of a broader attempt to describe systematically the effects of correlated
density fluctuations in dense polymer systems, both for static [12, 13, 44, 45] and dynam-
ical [29, 35, 46] properties. An important unresolved question is for instance whether the
predicted long-range repulsive forces of van der Waals type (“Anti-Casimir effect”) [13, 45]
are observable, for instance in the oscillatory decay of the standard density pair-correlation
function of dense polymer solutions. Since the results presented here challenge an important
concept of polymer physics, they should hopefully be useful for a broad range of theoreti-
cal approaches which commonly assume the validity of the Gaussian chain model down to
molecular scales [47, 48, 49]. This study shows that a polymer in dense solutions should not
be viewed as one soft sphere (or ellipsoid) [50, 51, 52], but as a hierarchy of nested segmental
correlation holes of all sizes aligned and correlated along the chain backbone (Fig. 2 (b)). We
note that similar deviations from Flory’s hypothesis have been reported recently for linear
polymers [16, 17, 47] and polymer gels and networks [53, 54]. The repulsive interactions
should also influence the polymer dynamics, since strong deviations from Gaussianity are
expected on the scale where entanglements become important, hence, quantitative predic-
tions for the entanglement length Ne have to be regarded with more care. The demonstrated
swelling of chains should be included in the popular primitive path analysis for obtaining Ne
[55], especially if ‘short’ chains (N < 500) are considered. The effect could be responsible
for observed deviations from Rouse behavior [26, 56] as may be seen by considering the cor-
relation function Cpq ≡ 〈Xp ·Xq〉 of the Rouse modes Xp = 1N
dnrn cos(npπ/N) where
p, q = 0, . . . , N − 1 [6, 57]. Using (rn − rm)2 = r2n + r2m − 2rn · rm for the segment size, this
correlation function can be readily expressed as an integral over the second moment of the
segmental size distribution
Cpq = −
(rn − rm)2
cos(npπ/N) cos(mpπ/N) (31)
which can be solved using our result Eq. (5). This implies for instance for p = q that
Cpp =
2(πp)2
1− π√
. (32)
The bracket entails an important correction with respect to the classical description given by
the prefactor [6]. We are currently working out how static corrections, such as those for Cpp,
may influence the dynamics for polymer chains without topological constraints. (This may
be realized, e.g., within the BFM algorithm by using the L26 moves described in Sec. IIIA.)
Moreover, for thin polymer films of width H the repulsive interactions are known to be
stronger than in the bulk [12]. This provides a mechanism to rationalize the trend towards
swelling observed experimentally [58] and confirmed computationally [21]:
= log(s)/H. (33)
(Prefactors omitted for clarity.) Here Rx(s) and bx denote the components of the segment
size and the effective bond length parallel to the film. It also explains the (at first sight
surprising) systematic increase of the polymer dynamics with decreasing film thickness [22].
Specifically, the parallel component of the monomer mean-squared displacement gx(t) is
expected to scale as gx(t) ≈ R2x(s(t)) ∝ t1/4(1 + log(t)/H) for long reptating chains where
s(t) ∝ t1/4 [6]. (The corresponding effect for the three-dimensional bulk should be small,
however.) For the same reason (flexible) polymer chains close to container walls must be
more swollen and, hence, faster on intermediate time scales than their peers in the bulk.
Acknowledgments
We thank T. Kreer, S. Peter and A.N. Semenov (all ICS, Strasbourg, France),
S.P. Obukhov (Gainesville, Florida) and M. Müller (Göttingen, Germany) for helpful discus-
sions. A generous grant of computer time by the IDRIS (Orsay) is also gratefully acknowl-
edged. J.B. acknowledges financial support by the IUF and from the European Community’s
“Marie-Curie Actions” under contract MRTN-CT-2004-504052.
APPENDIX A: MOMENTS OF THE SEGMENTAL SIZE DISTRIBUTION AND
THEIR GENERATING FUNCTION
Higher moments of the segmental size distribution G(r, s) can be systematically obtained
from its Fourier transformation
G(q, s) =
d3r G(r, s) exp(iq · r),
which is in this context sometimes called the “generating function” [59]. For an ideal
Gaussian chain, the generating function is then G0(q, s) = exp(−sq2a2) where we have
used a2 = b2/6 instead of the bond length b2 to simplify the notation. Moments of the
size distribution are given by proper derivatives of G(q, s) taken at q = 0. For example,
〈r2p〉 = (−1)p∆pG(q, s)q=0 (with ∆ being the Laplace operator with respect to the wave vec-
tor q). A moment of order 2p is, hence, linked to only one coefficient A2p in the systematic
expansion, G(q, s) =
p=0A2pq
2p, of G(q, s) around q = 0. For our example this implies
= (−1)p(2p+ 1)! A2p (A1)
in general and more specifically for a Gaussian distribution 〈r2p〉0 =
(2p+1)!
spa2p. The non-
Gaussian parameters read, hence,
αp(s) ≡ 1−
(2p+ 1)!
〈r2p〉
〈r2〉p
= 1− p! A2p
, (A2)
which implies (by construction) αp = 0 for a Gaussian distribution. As various moments of
the same global order 2p are linked to the same A2p they differ by a multiplicative constant
independent of the details of the (isotropic) distribution G(q, s). For example, 〈r2〉 = 6|A2|,
〈r4〉 = 120A4, 〈x2〉 = 〈y2〉 = 2|A2|, 〈x2y2〉 = 8A4 with x and y denoting the spatial
components of the segment vector r. Using Eq. (A2) for p = 2 it follows that
Kxy(s) ≡ 1−
〈x2y2〉
〈x2〉〈y2〉
= 1− 2A4
= α2(s), (A3)
i.e. the properties α2(s) and Kxy(s) discussed in Figs. 9 and 10 must be identical in general
provided that G(q, s) is isotropic and can be expanded in q2.
We turn now to specific properties ofG(q, s) computed for formally infinite polymer chains
in the melt. In practice, these results are also relevant for small segments in large chains,
N ≫ s ≫ 1, and, especially, for segments located far from the chain ends. These chains are
nearly Gaussian and the generating function can be written as G(q, s) = G0(q, s) + δG(q, s)
where δG(q, s) = −〈UG〉0 + 〈U〉0〈G〉0 is a small perturbation under the effective interaction
potential ṽ(q) given by Eq. (9). To compute the different integrals it is more convenient to
work in Fourier-Laplace space (q, t) with t being the Laplace variable conjugate to s:
δG(q, t) =
ds δG(q, s)e−st.
As illustrated in Fig. 14, there a three contributions to this perturbation: one due to in-
teractions between two monomers inside the segment (left panel), one due to interactions
between an internal monomer and an external one (middle panel) and one due to interac-
tions between two external monomers located on opposite sides (right panel). In analogy to
the derivation of the form factor described in Ref. [14] this yields:
δG(q, t) = −
(q2a2 + t)2
q2a2 + t−
(q2a2 + t)2
4πqa2ξ2
Arctan
a/ξ +
a/ξ +
q2a2 + t
q2a2 + t
4πqa4
Arctan
a/ξ +
a/ξ +
q2a2 + t
4πqa6
Arctan
q2a2 + t
. (A4)
The graph given in the left panel of Fig. 14 corresponds to the first two lines, the middle
panel to the third line and the right panel to the last one. Seeking for the moments we
expand δG(q, t) around q = 0. Having in mind chain strands counting many monomers
(s ≫ 1), we need only to retain the most singular terms for t → 0. Defining the two
dimensionless constants d = vξ/3πa4 = 12vξ/πb4 and c = (3π3/2a3ρ)−1 =
24/π3/b3ρ this
expansion can be written as
δG(q, t) = −
d a2q2 +
Γ(3/2)
c a2q2 + . . . (A5)
d a4q4 − 1
Γ(5/2)
c a4q4 + . . .
d a6q6 +
Γ(7/2)
c a6q6 + . . .
+ . . .
where we have used Euler’s Gamma function Γ(α) [60]. The first leading term at each order
in q2 — being proportional to the coefficient d — ensures the renormalization of the effective
bond length. The next term scaling with the coefficient c corresponds to the leading finite
strand size correction. Performing the inverse Laplace transformation Γ(α)/tα → sα−1 and
adding the Gaussian reference distribution G0(q, s) this yields the A2p-coefficients for the
expansion of G(q, s) around q = 0:
A0 = 1
A2 = −a2s
1 + d− c√
1 + 2d−
A6 = −
1 + 3d− 216
A8 = . . . (A6)
More generally, one finds
A2p =
(−1)p
(sa2)p
1 + pd− 3(2
pp!p)2
2(2p+ 1)!
From this result and using Eq. (A1) one immediately verifies that the moments of the
distribution are given by the Eqs. (11) and (12). Using Eq. (A2) one justifies similarly
Eq. (26) for the non-Gaussian parameter αp.
These moments completely determine the segmental distribution G(r, s) which is indi-
cated in Eq. (13). While at least in principle this may be done directly by inverse Fourier-
Laplace transformation of the correction δG(q, t) to the generating function it is helpful to
simplify further Eq. (A4). We observe first that δG(q, t) does diverge for strictly incom-
pressible systems (v → ∞) and one must keep v finite in the effective potential whenever
necessary to ensure convergence (actually everywhere but in the diagram corresponding to
the interaction between two external monomers). Since we are not interested in the wave
vectors larger than 1/ξ we expand δG(q, t) for ξ → 0 which leads to the much simpler
expression
δG(q, t) ≈ − vξq
3πa2(a2q2 + t)2
t(3a2q2 + t)
(a2q2 + t)2
Arctan[ aq√
+ o(vξ3). (A8)
The first term diverges as
v for diverging v. It renormalizes the effective bond length in
the zero order term which is indicated in the first line of Eq. (13). The next two terms scale
both as v0. Subsequent terms must all vanish for diverging v and can be discarded. It is
then easy to perform an inverse Fourier-Laplace transformation of the two relevant v0 terms.
This yields
δG(x, s) = G0(x, s)
with x = r/a
6n. This is consistent with the expression given in the second line of
Eq. (13).
We note finally that the intramolecular form factor F (q) = 1
n,m=1 〈exp(iq · (rn − rm)〉
of asymptotically long chains can be readily obtained from Eq. (A8). Observing that
〈exp(iq · (rn − rm)〉 =
d3r exp(iq · r)G(r, s) = G(q, s) one finds
δF (q) = 2
ds δG(q, s) = 2 δG(q, t = 0) = −2 vξ
, (A10)
where we used the third term of Eq. (A8) in the last step. The first term in Eq. (A8) is
discarded as before, since it renormalizes the effective bond length in the reference form
factor: F0(q) = 12/b
2q2 ⇒ 12/b2eq2. It follows, hence, that within first-order perturbation
theory
F (q) = F0(q) + δF (q) ≈ F0(q)
(A11)
as indicated by Eq. (30) in the Conclusion. This is equivalent to the result 1/F (q)−1/F0(q) ≈
q3/32ρ discussed in Refs. [14, 19] for polymer melts and anticipated by Schäfer [11] by
renormalization group calculations of semidilute solutions.
N nch τe Re Rg be(N)
6bg(N)
16 216 1214 11.7 4.8 2.998 2.939
32 215 3485 17.1 7.0 3.066 3.030
64 214 1.1 · 104 24.8 10.1 3.116 3.094
128 8192 3.3 · 104 35.6 14.5 3.153 3.139
256 4096 1.0 · 105 50.8 20.7 3.179 3.171
512 2048 3.2 · 105 72.2 29.5 3.200 3.193
1024 1024 1.0 · 106 103 42.0 3.216 3.212
2048 512 3.2 · 106 146 59.5 3.227 3.223
4096 256 9.7 · 106 207 85.0 3.235 3.253
8192 128 2.9 · 107 294 120 3.249 3.248
TABLE I: Various static properties of dense BFM melts of number density ρ = 0.5/8: the chain
length N , the number of chains nch per box, the relaxation time τe characterized by the diffusion
of the monomers over the end-to-end distance and corresponding to the circles indicated in Fig. 3,
the root-mean-squared chain end-to-end distance Re and the radius of gyration Rg of the total
chain (s = N − 1). The last two columns give estimates for the effective bond length from the
end-to-end distance, be(N) ≡ Re/(N − 1)1/2, and the radius of gyration, bg(N) ≡ Rg/
N . The
dashed line in Fig. 4 indicates be(N)
2. Apparently, both estimates increase monotonicly with N
reaching be(N) ≈
6bg(N) ≈ 3.2 for the largest chains available. Note that
6bg(N) < be(N) for
smaller N .
Property BFM BSM
Length unit lattice constant bead diameter
Temperature kBT 1 1
Number density ρ 0.5/8 0.84
Linear box size L 256 ≤ 62
Number of monomers nmon 1048576 ≤ 196608
Largest chain length N 8192 1024
Mean bond length 〈|ln|〉 2.604 0.97
l = 〈l2n〉1/2 2.636 0.97
Effective bond length be 3.244 1.34
ρb3e 2.13 2.02
C∞ = (be/l)
2 1.52 1.91
lp = l(C∞ + 1)/2 3.32 1.41
24/π3/ρb3e 0.41 0.44
c1/ce 1.0 1.2
c2/ce 1.0 1.1
c3/ce 1.0 1.0
c4/ce 1.1 1.2
c5/ce 1.1 0.9
cP = c1(be/l)
2/8 0.078 0.124
Dimensionless compressibility g 0.245 0.08
Compression modulus v ≡ 1/gρ 66.7 14.9
vρ/b3eρ 0.96 1.8
TABLE II: Comparison of some static properties of dense BFM and BSM melts. The first six rows
indicate conventions and operational parameters. The effective bond length be and the swelling
coefficients cp (defined in Eq. (5)) are determined from the first five even moments of the segmental
size distribution. The dimensionless compressibility g = S(q → 0)/ρ has been obtained from the
total static structure factor S(q) = 1
∑nmon
k,l=1 〈exp(iq · (rk − rl))〉 in the zero wave vector limit
as shown at the end of Ref. [14]. The values indicated correspond to the asymptotic long chain
behavior. Properties of very small chains deviate slightly.
I ~ −3 I ~ −9
I ~ 45
I ~ −9−
r   n m
FIG. 1: (Color online) Sketch of a polymer chain of length N in a dense melt in d = 3 dimensions.
As notations we use ri for the position vector of a monomer i, li = ri+1 − ri for its bond vector,
r = rm−rn for the end-to-end vector of the chain segment between the monomers n and m = n+s
and r = ||r|| for its length. Segment properties, such as the 2p-th moments
, are averaged over
all possible pairs of monomers (n,m) of a chain and over all chains. The second moment (p = 1)
is denoted Re(s) =
, the total chain end-to-end distance is Re(s = N − 1). The dashed
lines show the relevant graphs of the analytical perturbation calculation outlined in Sec. II B.
The numerical factors indicate for infinite chains (without chain end effects) the relative weights
contributing to the 1/
s-swelling of Re(s) indicated in Eq. (10).
a bdensity ρ
c(r,s)
correlation hole
Segmental
= const
R(s) Repulsion
FIG. 2: (Color online) Role of incompressibility and chain connectivity in dense polymer solutions
and melts. (a) Sketch of the segmental correlation hole of a marked chain segment of curvilin-
ear length s. Density fluctuations of chain segments must be correlated, since the total density
fluctuations (dashed line) are small. Consequently, a second chain segment feels an entropic re-
pulsion when both correlation holes start to overlap. (b) Self-similar pattern of nested segmental
correlation holes of decreasing strength u(s) ≈ s/ρR(s)3 ≈ ce/
s aligned along the backbone of
a reference chain. The large dashed circle represents the classical correlation hole of the total
chain (s ≈ N) [5]. This is the input of recent approaches to model polymer chains as soft spheres
[50, 52]. We argue that incompressibility on all scales and chain connectivity leads to a short
distance repulsion of the segmental correlation holes, which increases with decreasing s.
L06 (conserved topology)
L26,SS
L26,SS,DB
FIG. 3: (Color online) Diffusion time τe over the (root-mean-squared) chain end-to-end distance
Re(N − 1) as a function of chain length N for different versions of the Bond Fluctuation Model
(BFM). All data indicated are for the high number density (ρ = 0.5/8) corresponding to a polymer
melt with half the lattice sites being occupied. We have obtained τe = R
e(N−1)/6Ds from the self-
diffusion coefficient Ds measured from the free diffusion limit of the mean-squared displacement
of all monomers
δr(t)2
= 6Dst. Data from the classical BFM with topology conserving local
Monte Carlo (MC) moves in 6 spatial directions (L06) [26] are represented as stars. All other data
sets use topology violating local MC moves in 26 lattice directions (L26). If only local moves are
used, L26-dynamics is even at relatively short times perfectly Rouse like which allows the accurate
determination of Ds although the monomers possibly have not yet moved over Re(N − 1) for the
largest chain lengths considered. Additional slithering snake (SS) moves increase the efficiency of
the algorithm by approximately an order of magnitude (squares,bold line). The power law exponent
is changed from 2 to an empirical 1.62 (dashed line) if in addition we perform double-bridging (DB)
moves.
N=128
N=256
N=512
N=1024
N=2048
N=4096
N=8192
=3.244
Eq.(5)
FIG. 4: (Color online) Mean-squared segment size Re(s)
2/s vs. curvilinear distance s. We present
BFM data for different chain length N at number density ρ = 0.5/8. The averages are taken over
all possible monomer pairs (n,m = n+s). The statistics deteriorates, hence, for large s. Log-linear
coordinates are used to emphasize the power law swelling over several orders of magnitude of s.
The data approach the asymptotic limit (horizontal line) from below, i.e. the chains are swollen.
This behavior is well fitted by Eq. (5) for 1 ≪ s ≪ N (bold line). Non-monotonic behavior is found
for s → N , especially for small N . The dashed line indicates the measured total chain end-to-end
distances, be(N)
2 ≡ Re(N − 1)2/(N − 1) from Tab. I, showing even more pronounced deviations
from the asymptotic limit. The dash-dotted line compares this data with Eq. (19).
=3.235
=3.240
=3.244
=3.250
=3.255
BFM ρ=0.5/8, N=2048: 
=3.244, c
=0.41=c
FIG. 5: (Color online) Replot of the mean-squared segment size as y = K1(s) = 1 − Re(s)2/b2es
vs. x = c1/
s, as suggested by Eq. (5), for different trial effective bond lengths be as indicated.
Only BFM chains of length N = 2048 are considered for clarity. This procedure is very sensitive
to the value chosen and allows for a precise determination. It assumes, however, that higher order
terms in the expansion of K1(s) may be neglected. The value be is confirmed from a similar test
for higher moments (Fig. 6).
x ~1/s
BFM b
=3.244
BSM b
=1.34
too small s !
FIG. 6: (Color online) Critial test of Eq. (5) where the rescaled moments y = Kp(s) of the segment
size distribution (defined in Eq. (1)) are plotted vs. x =
3(2pp!p)2
2(2p+1)!
. We consider the first five even
moments (p = 1, . . . , 5) for the BFM with N = 2048 and the BSM with N = 1024. Also indicated
is the rescaled radius of gyration, y = 5/8 (1 − 6R2g(s)/b2e(s + 1)), as a function of x = c1/
(filled circles). The BSM data has been shifted upwards for clarity. Without this shift a perfect
data collapse is found for both models and all moments. Keeping the same effective bond length
be for all moments of each model we fit for the swelling coefficients cp by rescaling the horizontal
axis. We find be ≈ 3.244 for the BFM and 1.34 for the BSM. If be is chosen correctly, all data
sets extrapolate linearly to zero for large s (x → 0). The swelling coefficients found are close the
theoretical prediction ce, as indicated in Tab. II. The plot demonstrates that the non-Gaussian
deviations scale as the segmental correlation hole, u(s) ∼ ce/
s and this for all moments as long
as x ≪ 1. The saturation at large x is due to the finite extensibility of short chain segments. Since
this effect becomes more marked for larger moments, the fit of be is best performed for p = 1.
u(s)  c
Kλ(s)
λ=2: N=1024
λ=2: N=2048
λ=2: N=4096
λ=4: N=2048
λ=8: N=2048
λ=16: N=2048
Insufficient statistics & Finite chain length effects
⇒ ω = 3/2
FIG. 7: (Color online) Plot of Kλ(s) as a function of u(s)c1/ce ∼ 1/
s using the measured
u(s) ≡
24/π3s/ρRe(s)
3. For λ = 2 (corresponding to two segments being connected) BFM and
BSM data are compared. Several λ values are given for N = 2048 BFM chains. For chain segments
with 1 ≪ s ≪ N all data sets collapse on the bisection line confirming the so-called “recursion
relation” Kλ ≈ u proposed by Semenov and Johner [12]. The statistics becomes insufficient for
large s (left bottom corner). Systematic deviations arise for s → N due to additional finite-N
effects.
N=128
N=512
N=2048
N=4096
N=8192
Slope ω=3/2
Chain end effects
FIG. 8: (Color online) The bond-bond correlation function P (s)/cP as a function of the curvilinear
distance s. Various chain lengths are given for BFM. Provided that 1 ≪ s ≪ N , all data sets
collapse on the power law slope with exponent ω = 3/2 (bold line) as predicted by Eq. (23). The
dash-dotted curve P (s) ≈ exp(−s/1.5) shows that exponential behavior is only compatible with
very small chain lengths. The dashed lines correspond to the theoretical prediction, Eq. (24), for
short chains with N = 16, 32, 64 and 128 (from left to right).
BFM N=1024
BFM N=2048
BFM N=4096
BFM N=8192
BSM N=1024
y = 6x/5
y = 111x/35
y = 605x/105
small s !
Noise !
FIG. 9: (Color online) Non-Gaussian parameter αp(s) computed for the end-to-end distance of
chain segments as a function of ce/
s. Perfect data collapse for all chain lengths and both sim-
ulation models is obtained for each p. A linear relationship over nearly two orders of magnitude
is found as theoretically expected. Data for three moments (p = 2, 3, 4) are indicated showing
a systematic increase of non-Gaussianity with p. The data curvature for small s becomes more
pronounced for larger p.
N=256
N=512
N=1024
N=2048
N=4096
N=8192
m=n+s
FIG. 10: (Color online) Plot of Kxy(s) = 1−
averaged over all pairs of monomers
(n,m = n+s) and three different direction pairs as a function of ce/
s. As indicated by the sketch
at the bottom of the figure, Kxy(s) measures the correlation of the components of the segment
vector r. All data points collapse and show again a linear relationship Kxy ≈ u(s). Different
directions are therefore coupled! No curvature is observed over two orders of magnitude confirming
that higher order perturbation corrections are negligible. Noise cannot be neglected for large
s > 100 and finite segment-size effects are visible for s ≈ 1.
0.0 0.5 1.0 1.5 2.0
n=r/b
BFM s=32
BFM s=64
BFM s=128
BSM s=32
BSM s=64
BSM s=128
GaussEnhancem
entDepletion
FIG. 11: (Color online) Segment size distribution y = G(r, s)(bes
1/2)3 vs. n = r/bes
1/2 for several
s as indicated in the figure. Only data for BFM with N = 2048 and BSM with N = 1024 are
presented. (A similar plot can be achieved by renormalizing the axes using Re(s) instead of bes
1/2).
The bold line denotes the Gaussian behaviour y = (3/2π)3/2 exp(−3n2/2). One sees that compared
to this reference the measured distributions are depleted for small n ≪ 1 (where the data does not
scale) and enhanced for n ≈ 1.
0.0 0.5 1.0 1.5 2.0 2.5
n=r/b
BFM s=31
BFM s=63
BFM s=127
BSM s=32
BSM s=64
BSM s=128
f(n) from Eq. (3)
Enhancement
epletion
FIG. 12: (Color online) Deviation δG(r, s) = G(r, s) − G0(r, s) of the measured segmental
size distribution from the Gaussian behavior G0(r, s) expected from Flory’s hypothesis for sev-
eral s and both models as indicated in the figure. As suggested by Eq. (3), we have plotted
y = (δG(r, s)/G0(r, s))/(ce/
s) as a function of n = r/be
s. The Gaussian reference distribution
has been computed according to Eq. (2) for the measured effective bond length be. A close to
perfect data collapse is found for both models. This shows that the deviation scales linearly with
u(s) ≈ ce/
s, as expected. The bold line indicates the universal function of f(n) predicted by
Eq. (3).
0.0 0.5 1.0 1.5 2.0 2.5
m=r/R
BFM s=31
BFM s=63
BFM s=127
BSM s=32
BSM s=64
BSM s=128
f(m) from Eq. (29)
Enhancement
letio
FIG. 13: (Color online) Replot of the relative deviation of the measured segment size distribution,
y = (δG(r, s)/G0(r, s)))/(ce/
s), as a function of m = r/Re(s). The figure highlights that the
measured segment size is the only length scale relevant for describing the deviation from Flory’s
hypothesis. The same data sets and symbols are used as in the previous Fig. 12.
v(k) v(k) v(k)~ ~ ~
q q q q q
q−k q−k
FIG. 14: Interaction diagrams used in reciprocal space for the calculation of δG(q, t) in the scale
free limit. There exist three nonzero contributions to first-order perturbation, the first involving
two points inside the segment (first two lines of Eq. (A4)), the second one point inside and one
outside the segment (third line of Eq. (A4)) and the third one point on either side of the segment
(last line of Eq. (A4)). Momentum q flows from one correlated point to the other. Integrals are
performed over the momentum k. Dotted lines denote the effective interactions ṽ(k) given by
Eq. (9), bold lines the propagators which carry each a momentum q or q − k as indicated.
[1] M. Rubinstein and R. Colby, Polymer Physics (Oxford University Press, 2003).
[2] P. J. Flory, J. Chem. Phys. 13, 453 (1945).
[3] P. J. Flory, J. Chem. Phys. 17, 303 (1949).
[4] P. J. Flory, Statistical Mechanics of Chain Molecules (Oxford University Press, New York,
1988).
[5] P. G. de Gennes, Scaling Concepts in Polymer Physics (Cornell University Press, Ithaca, New
York, 1979).
[6] M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Clarendon Press, 1986).
[7] S. Edwards, Proc. Phys. Soc. 85, 613 (1965).
[8] S. Edwards, Proc. Phys. Soc. 88, 265 (1966).
[9] S. Edwards, J.Phys.A: Math.Gen. 8, 1670 (1975).
[10] M. Muthukumar and S. Edwards, J. Chem. Phys. 76, 2720 (1982).
[11] L. Schäfer, Excluded Volume Effects in Polymer Solutions (Springer-Verlag, New York, 1999).
[12] A. N. Semenov and A. Johner, Eur. Phys. J. E 12, 469 (2003).
[13] A. N. Semenov and S. P. Obukhov, J. Phys.: Condens. Matter 17, 1747 (2005).
[14] P. Beckrich, A. Johner, A. N. Semenov, S. P. Obukhov, H. C. Benôıt, and J. P. Wittmer,
Macromolecules (2007), in press, cond-mat/0701261.
[15] P. Beckrich, Ph.D. thesis, Université Louis Pasteur, Strasbourg, France (2006).
[16] L. Schäfer, M. Müller, and K. Binder, Macromolecules 33, 4568 (2000).
[17] R. Auhl, R. Everaers, G. Grest, K. Kremer, and S. Plimpton, J. Chem. Phys. 119, 12718
(2003).
[18] J. P. Wittmer, H. Meyer, J. Baschnagel, A. Johner, S. P. Obukhov, L. Mattioni, M. Müller,
and A. N. Semenov, Phys. Rev. Lett. 93, 147801 (2004), cond-mat/0404457.
[19] J. P. Wittmer, P. Beckrich, A. Johner, A. N. Semenov, S. P. Obukhov, H. Meyer, and
J. Baschnagel, Europhys. Lett. 77, 56003 (2007), cond-mat/0611322.
[20] J. P. Wittmer, P. Beckrich, F. Crevel, C. C. Huang, A. Cavallo, T. Kreer, and H. Meyer,
Computer Physics Communications (2007), in press, cond-mat/0610359.
[21] A. Cavallo, M. Müller, J. P. Wittmer, and A. Johner, J. Phys.: Condens. Matter 17, 1697
(2005), cond-mat/0412373.
[22] H. Meyer, T. Kreer, A. Cavallo, J. P. Wittmer, and J. Baschnagel, Eur. Phys. J. Special Topics
141, 167 (2007), cond-mat/0609127.
[23] J. Baschnagel, J. P. Wittmer, and H. Meyer, in Computational Soft Matter: From Synthetic
Polymers to Proteins, edited by N. Attig (NIC Series, Jülich, 2004), vol. 23, pp. 83–140,
cond-mat/0407717.
[24] I. Carmesin and K. Kremer, Macromolecules 21, 2819 (1988).
[25] H. Deutsch and K. Binder, J. Chem. Phys 94, 2294 (1991).
[26] W. Paul, K. Binder, D. Heermann, and K. Kremer, J. Phys. II 1, 37 (1991).
[27] A. Kron, Polym. Sci. USSR 7, 1361 (1965).
[28] F. Wall and F. Mandel, J. Chem. Phys. 63, 4592 (1975).
[29] L. Mattioni, J. P. Wittmer, J. Baschnagel, J.-L. Barrat, and E. Luijten, Eur. Phys. J. E 10,
369 (2003), cond-mat/0212433.
[30] N. Karayiannis, A. Giannousaki, V. Mavrantzas, and D. Theodorou, J. Chem. Phys. 117,
5465 (2002).
[31] B. J. Banaszak and J. J. de Pablo, J. Chem. Phys. 119, 2456 (2003).
[32] M. Müller, J. P. Wittmer, and J.-L. Barrat, Europhys. Lett. 52, 406 (2000), cond-
mat/0006464.
[33] T. Kreer, J. Baschnagel, M. Müller, and K. Binder, Macromolecules 34, 1105 (2001).
[34] J. M. Deutsch, Phys. Rev. Lett. 54, 56 (1985).
[35] A. N. Semenov, in Theoretical Challenges in the Dynamics of Complex Fluids, edited by
T. McLeish (Kluwer, Dordrecht, 1997), pp. 63–86.
[36] H. Meyer and F. Müller-Plathe, J. Chem. Phys. 115, 7807 (2001).
[37] H. Meyer and F. Müller-Plathe, Macromolecules 35, 1241 (2002).
[38] K. Kremer and G. Grest, J. Chem. Phys. 92, 5057 (1990).
[39] M. Allen and D. Tildesley, Computer Simulation of Liquids (Oxford University Press, Oxford,
1994).
[40] O. Kratky and G. Porod, Rec. Trav. Chim. 68, 1106 (1949).
[41] J. Hansen and I. McDonald, Theory of simple liquids (Academic Press, New York, 1986).
[42] B. Alder and T. Wainwright, Phys. Rev. A 1, 18 (1970).
[43] P. Beckrich, J. P. Wittmer, H. Meyer, S. P. Obukhov, A. N. Semenov, and A. Johner (2007),
in preparation.
[44] A. N. Semenov, Journal de Physique II France 6, 1759 (1996).
[45] S. P. Obukhov and A. N. Semenov, Phys. Rev. Lett. 95, 038305 (2005).
[46] A. N. Semenov and M. Rubinstein, Eur. Phys. J. B 1, 87 (1998).
[47] J. Curro, K. S. Schweizer, G. S. Grest, and K. Kremer, J. Chem. Phys. 91, 1359 (1991).
[48] M. Fuchs, Z. Phys. B 5, 521 (1997).
[49] M. Schulz, H. Frisch, and P. Reineker, New Journal of Physics 6, 77/1 (2004).
[50] C. Likos, Physics Reports 348, 267 (2001).
[51] F. Eurich and P. Maass, J. Chem. Phys. 114, 7655 (2001).
[52] G. Yatsenko, E. J. Sambriski, M. A. Nemirovskaya, and M. Guenza, Phys. Rev. Lett. 93,
257803 (2004).
[53] J.-U. Sommer and K. Saalwächter, Eur. Phys. J. E 18, 167 (2005).
[54] C. Svaneborg, G. Grest, and R. Everaers, Europhys. Lett. 72, 760 (2005).
[55] R. Everaers, S. Sukumaran, G. Grest, C. Svaneborg, A. Sivasubramanian, and K. Kremer,
Science 303, 823 (2004).
[56] W. Paul and G. Smith, Rep. Progr. Phys. 67, 1117 (2004).
[57] P. Verdier, J. Chem. Phys. 45, 2118 (1966).
[58] R. Jones, S. Kumar, D. Ho, R. Briber, and T. Russell, Nature 400, 146 (1999).
[59] N. G. van Kampen, Stochastic processes in physics and chemistry (North-Holland, Amster-
dam, 1992).
[60] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions (Dover, New York,
1964).
[61] It is interesting to compare the numerical value I(∞) ≈ 1.59 obtained for the r.h.s of Eq. (16)
with the coefficients one would obtain by computing Eq. (15) either with the effective potential
ṽ(q) for infinite chains given by Eq. (9) or with the Padé approximation, Eq. (17). Within
these approximations of the full linear response formula, Eq. (14), the coefficients can be
obtained directly without numerical integration yielding overall similar values. In the first
case we obtain 15/8 ≈ 1.87 and in the second 11/8 ≈ 1.37. While the first value is clearly
not compatible with the measured end-to-end distances, the second yields a reasonable fit,
especially for small N < 1000, when the data is plotted as in Fig. 5. Ultimately, for very long
chains the correct coefficient should be 1.59 as is indicated by the dash-dotted line in Fig. 4.
[62] These topology non-conserving moves yield configurations which are not accessible with the
classical scheme with jumps in 6 directions only. Concerning the static properties we are
interested in this paper both system classes are practically equivalent. This has been confirmed
by comparing various static properties and by counting the number of monomers which become
“blocked” (in absolute space or with respect to an initial group of neighbor monomers) once
one returns to the original local scheme. Typically we find about 10 blocked monomers for a
system of 220 monomers. The relative difference of microstates is therefore tiny and irrelevant
for static properties. Care is needed, however, if the equilibrated configurations are used to
investigate the dynamics of the topology conserving BFM version. The same caveats arise for
the slithering snake and double-bridging moves.
[63] We have verified that the alternative definition 〈em=n+s · en〉 with en being the normalized
bond vector yields very similar results. This is due to the weak bond length fluctuations,
specifically at high densities, in both coarse-grained models under consideration. It is possi-
ble that other models show slightly different power law amplitudes cP depending on which
definition is taken.
	Flory's ideality hypothesis revisited
	Physical idea and sketch of the perturbation calculation
	Scaling arguments
	Perturbation calculation
	Computational models and technical details
	Bond fluctuation model
	Bead spring model
	Numerical results
	The swelling of chain segments
	Chain connectivity and recursion relation
	Intrachain bond-bond correlations 
	Higher moments and associated coefficients
	The segmental size distribution
	Conclusion 
	Acknowledgments
	Moments of the segmental size distribution and their generating function
	References

  Presenting theoretical arguments and numerical results we demonstrate
long-range intrachain correlations in concentrated solutions and melts of long
flexible polymers which cause a systematic swelling of short chain segments.
They can be traced back to the incompressibility of the melt leading to an
effective repulsion $u(s) \approx s/\rho R^3(s) \approx ce/\sqrt{s}$ when
connecting two segments together where $s$ denotes the curvilinear length of a
segment, $R(s)$ its typical size, $ce \approx 1/\rho be^3$ the ``swelling
coefficient", $be$ the effective bond length and $\rho$ the monomer density.
The relative deviation of the segmental size distribution from the ideal
Gaussian chain behavior is found to be proportional to $u(s)$. The analysis of
different moments of this distribution allows for a precise determination of
the effective bond length $be$ and the swelling coefficient $ce$ of
asymptotically long chains. At striking variance to the short-range decay
suggested by Flory's ideality hypothesis the bond-bond correlation function of
two bonds separated by $s$ monomers along the chain is found to decay
algebraically as $1/s^{3/2}$. Effects of finite chain length are considered
briefly.

Introduction, one sees that both equations are essentially identical —
taken apart, however, that they depend on b and be. Note the conspicuous factor (b/be)
in Eq. (12) which would strongly reduce the empirical swelling coefficients cp = ce(b/be)
for large p if b and be were different.
Interpretation of first-loop results in different contexts. The above perturbation results
may be used directly to describe the effect of a weak excluded volume v on a reference system
of perfectly ideal polymer melts with Kuhn segment length b where all interactions have been
switched off (v = 0). It is expected to give a good estimation for the effective bond length
be only for a small Ginzburg parameter: Gz ≪ 1. For the dense melts we want to describe
this does not hold (Sec. III) and one cannot hope to find a good quantitative agreement
with Eq. (11). Note also that large wave vectors contribute strongly to the leading Gaussian
term. The effective bond length be is, hence, strongly influenced by local and non-universal
effects and is very difficult to predict in general.
Our much more modest goal is to predict the coefficient of the 1/
s-perturbation and
to express it in terms of a suitable variational reference Hamiltonian characterized by a
conveniently chosen Kuhn segment b and the measured effective bond length be (instead of
Eq. (11)). Following Muthukumar and Edwards [10], we argue that for dense melts b should
be renormalized to be to take into account higher order graphs. No strict mathematical
proof can be given at present that the infinite number of possible graphs must add up in
this manner. Our hypothesis relies on three observations:
• The general scaling argument discussed in Sec. IIA states that we have only one
relevant length scale in this problem, the typical segment size R(s) ≈ be
s itself. The
incompressibility constraint cannot generate an additional scale. It is this size R(s)
which sets the strength of the effective interaction which then in turn feeds back to
the deviations of R(s) from Gaussianity. Having a bond length b in addition to the
effective bond length be associated with R(s) would imply incorrectly a second length
scale b
s varying independently with the bulk modulus v. (We will check explicitly
below in Fig. 13 that there is only one length scale.) This implies b/be = const v
• Thus, since by construction b/be = 1 for v → 0, it follows that both lengths should be
equal for all v.
• We know from Eq. (12) that the empirical coefficients cp = ce(b/be)2p−3 should depend
strongly on the moment considered if the ratio b/be is not close to unity. It will be
shown below (Fig. 6) that cp/ce ≈ 1 for all p. This implies b ≈ be.
Finite chain size effects. To describe properly finite chain size corrections Eq. (9) must
be replaced by the general linear response formula
ṽ(q)
+ ρF (q) (14)
with F (q) = NfD(x) being the form factor of the Gaussian reference chain given by Debye’s
function fD(x) = 2(e
−x − 1 + x)/x2 with x = (qb)2N/6 [6]. This approximation allows
in principle to compute, for instance, the (mean-squared) total chain end-to-end distance,
A = (rN −r1)2. One verifies readily (see [6], Eq. (5.III.9)) that the effect of the perturbation
may be expressed as
〈A〉0 〈U〉0 − 〈AU〉0 =
(2π)3
ṽ(q)
ds s2(N − s) exp
q2b2s
. (15)
We take now first the integral over s. In the remaining integral over q small q wave vectors
contribute to the
N -swelling while large q renormalize the effective bond length of the
dominant Gaussian behaviour linear in N (as discussed above). Since we wish to determine
the non-Gaussian corrections, we may focus on small wave vectors q ≪ 1/ξ. Since in this
limit 1/v = ρg ≪ ρF (q), one can neglect in Eq. (14) the 1/v contribution to the inverse
effective interaction potential. We thus continue the calculation using the much simpler
ṽ(N, x) = 1/(NρfD(x)). This allows us to express the swelling as
〈(rN − r1)2〉
I(xu). (16)
To simplify the notation we have set here finally b = be in agreement with the hypothesis
discussed above. The numerical integral I(xu) =
dx . . . over x is slowly convergent at
infinity. As a consequence the estimate I(∞) = 1.59 may be too large for moderate chain
lengths. In practice, convergence is not achieved for values xu(N) ≈ (b/ξ)2N corresponding
to the screening length ξ.
We remark finally that numerical integration can be avoided for various properties if the
Padé approximation of the form factor, F (q) = N/(1 + (qb)2N/12), is used. This allows
analytical calculations by means of the simplified effective interaction potential
ṽ(q) =
12ρb3
. (17)
This has been used for instance for the calculation of finite chain size effects for the bond-
bond correlation function discussed in Sec. IVC below [61].
III. COMPUTATIONAL MODELS AND TECHNICAL DETAILS
A. Bond fluctuation model
A widely-used lattice Monte Carlo scheme for coarse-grained polymers. The body of
our numerical data comes from the three dimensional bond fluctuation model (BFM) —
a lattice Monte Carlo (MC) algorithm where each monomer occupies eight sites of a unit
cell of a simple cubic lattice [24, 25, 26]. Our version of the BFM with 108 bond vectors
corresponds to flexible athermal chain configurations [23]. All length scales are given in units
of the lattice constant and time in units of Monte Carlo Steps (MCS). We use cubic periodic
simulation boxes of linear size L = 256 containing nmon = ρL
3 = 220 ≈ 106 monomers.
This monomer number corresponds to a monomer number density ρ = 0.5/8 where half of
the lattice sites are occupied (volume fraction 0.5). The large system sizes used allow us to
suppress finite box-size effects for systems with large chains. Using a mix of local, slithering
snake [27, 28, 29], and double-bridging [17, 23, 30, 31] MC moves we were able to equilibrate
dense systems with chain lengths up to N = 8192.
Equilibration and sampling of high-molecular BFM melts. Standard BFM implementa-
tions [26, 32, 33] use local MC jumps to the 6 closest lattice sites to prevent the crossing
of chains and conserve therefore the chain topology. These “L06” moves lead to very large
relaxation times, scaling at least as τe ∼ N3, as may be seen from Fig. 3 (stars). The relax-
ation time τe = R
e/6Ds indicated in this figure has been estimated from the self-diffusion
coefficient Ds obtained from the mean-square displacements of all monomers in the free
diffusion limit. (For the largest chain indicated for L06 dynamics only a lower bound for
τe is given.) Instead of this more realistic but very slow dynamical scheme we make jump
attempts to the 26 sites of the cube surrounding the current monomer position (called “L26”
moves). This allows the chains to cross each other which dramatically speeds up the dynam-
ics, especially for long chains (N > 512). If only local moves are considered, the dynamics
is perfectly consistent with the Rouse model [6]. As shown in Fig. 3, we find τe ≈ 530N2 for
L26 dynamics. This is, however, still prohibitive by large for sampling configurations with
the longest chain length N we aim to characterize [62].
Slithering snake moves. In addition to the local moves one slithering snake move per
chain is attempted on average per MCS corresponding to the displacement of N monomer
along the chain backbone. Note that in our units two spatial displacement attempts per
MCS are performed on average per monomer, one for a local move and one for a snake
move. (In practice, it is computationally more efficient for large N to take off a monomer
at one chain end and to paste it at the other leaving all other monomers unaltered. Before
dynamical measurements are performed the original order of beads must then be restored.)
Interestingly, a significantly larger slithering snake attempt frequency would not be useful
since the relaxation time of slithering snakes without or only few local moves increases
exponentially with mass [29, 34] due to the correlated motion of snakes [35]. In order to
obtain an efficient free snake diffusion (with a chain length independent curvilinear diffusion
coefficient Dc(N) ∼ N0 and τe ≈ N2/Dc(N) ∼ N2 [28, 29]) it is important to relax density
fluctuations rapidly by local dynamical pathways. As shown in Fig. 3 (squares), we find
a much reduced relaxation time τe ≈ 40N2 which is, however, still unconveniently large
for our longest chains. Note that most of the CPU time is still used by local moves. The
computational load per MCS remains therefore essentially chain length independent.
Advantages and pitfalls of double-bridging moves. Double-bridging (DB) moves are very
useful for high densities and help us to extend the accessible molecular masses close to
104. As for slithering snake moves we use all 108 bond vectors to switch chain segments
between two different chains. Only chain segments of equal length are swapped to conserve
monodispersity. Topolocial constraints are again systematically and deliberately violated.
Since more than one swap partner is possible for a selected first monomer, delicate detailed
balance questions arise. This is particularly important for short chains and is discussed
in detail in Ref. [23]. Technically, the simplest solution to this problem is to refuse all
moves with more than one swap partner (to be checked both for forward and back move).
The configurations are screened with a frequency fDB for possible DB moves where we
scan in random order over the monomers. The frequency should not be too large to avoid
(more or less) immediate back swaps and monomers should move at least out of the local
monomer cage and over a couple of lattice sites. We use fDB = 0.1 between DB updates for
the configurations reported here. (The influence of fDB on the performance has not been
explored systematically, but preliminary results suggest a slightly smaller DB frequency for
future studies.) The diffusion times over the end-to-end distance for this case are indicated
in Tab. I. As shown in Fig. 3, we find empirically τe(N) ≈ 13N1.62. For N = 8192 this
corresponds to 3 · 107 MCS. This allows us even for the largest chain lengths to observe
monomer diffusion over several Re within the 10
8 MCS which are feasible on our XEON-PC
processor cluster.
The efficiency of DB moves is commonly characterized in terms of the relaxation time
τee of the end-to-end vector correlation function [30, 31]. For normal chain dynamics this
would indeed characterize the longest relaxation time of the system, i.e. τe ≈ τee. For the
double-bridging this is, however, not sufficient since density fluctuations do not couple to the
bridging moves and can not be relaxed. We find therefore that configurations equilibrate
on time scales given by τe rather than by τee ≪ τe. This may be verified, for instance,
from the time needed for the distribution Re(s) (and especially its spatial components) to
equilibrate. The criterion given in the literature [30] is clearly not satisfactory and may lead
to insufficiently equilibrated configurations. In summary, equilibration with DB moves still
requires monomer diffusion over the typical chain size, however at a much reduced price.
Some properties of our configurations. The Tables I and II summarize some system
properties obtained for our reference density ρ = 0.5/8. Averages are performed over all
chains and 1000 configurations. These configurations may be considered to be independent
for N < 4096. Only a few independent configurations exist for the largest chain length
N = 8192 which has to be considered with some care. Taking apart this system, chains are
always much smaller than the box size. For asymptotically long chains, we obtain an average
bond length 〈|l|〉 ≈ 2.604, a root-mean-squared bond length l ≡ 〈l2〉1/2 ≈ 2.635 and an
effective bond length be ≈ 3.244 — as we will determine below in Sec. IVA. This corresponds
to a ratio C∞ ≡ b2e/l2 ≈ 1.52 and, hence, to a persistence length lp = l(C∞ + 1)/2 ≈ 3.32
[23]. Especially, we find from the zero wave vector limit of the total structure factor S(q)
a low (dimensionless) compressibility g = S(q → 0)/ρ ≈ 0.246 which compares well with
real experimental melts. From the measured bulk compression modulus v ≡ 1/g(ρ)ρ and
the effective bond length be one may estimate a Ginzburg parameter Gz =
vρ/b3eρ ≈ 0.96.
Following Ref. [13] the interaction parameter v is supposed here to be given by the full
inverse compressibility and not just by the second virial coefficient.
B. Bead spring model
Hamiltonian. Additionally, molecular dynamics simulations of a bead-spring model
(BSM) [36] were performed to dispel concerns that our results are influenced by the un-
derlying lattice structure of the BFM. The model is derived from a coarse-grained model
for polyvinylalcohol which has been employed to study polymer crystallization [37]. It is
characterized by two potentials: a non-bonded potential of Lennard-Jones (LJ) type and a
harmonic bond potential. While the often employed Kremer-Grest model [38] uses a 12− 6
LJ potential to describe the non-bonded interactions Unb(r), our non-bonded potential has
a softer repulsive part. It is given by
Unb(r) = 1.511
, (18)
which is truncated and shifted at the minimum at rmin ≈ 1.15. Note that all length scales
are given in units of σ0 and we use LJ units [39] for all BSM data (mass m = 1, Boltzmann
constant kB = 1). The parameters of the bond potential, Ub(r) = 1120(r− lb)2, are adjusted
so that the average bond length l(ρ = 0.84) ≈ lb = 0.97 is approximately the same as in the
standard Kremer-Grest model [38]. The average bond length and the root-mean-squared
bond length are almost identical for the BSM due to the very stiff bond potential. Since
rmin/l ≈ 1.16 bonded monomers penetrate each other significantly.
Equilibration and sampling. We perform standard molecular dynamics simulations in the
canonical ensemble with a Langevin thermostat (friction constant Γ = 0.5) at temperature
T = 1. The equations of motion are integrated by the velocity-Verlet algorithm [39]. To
improve the statistics for large chain length, we have implemented additional double-bridging
moves. Since only few of these MC moves are accepted per unit time, this does affect neither
the stability nor the accuracy of the molecular dynamics sweeps.
Some properties obtained. For clarity, we show only data for chain length N = 1024
and number density ρ = 0.84, the typical melt density of the Kremer-Grest model [38].
For the reported data we use periodic simulation boxes of linear size L ≈ 62 containing
nmon = 196608 monomers, but we have also sampled different boxes sizes (up to L = 77.5)
to check for finite box-size effects. For the reference density a dimensionless compressibility
g ≈ 0.08 is found which is about three times smaller than for our BFM melt. For the
effective bond length we obtain be ≈ 1.34, i.e. BSM chains (C∞ ≈ 1.91, lp ≈ 1.41) are
slightly stiffer than the corresponding BFM polymers. Fortunately, the product ρb3e ≈ 2 is
roughly similar in both models and one expects from Eq. (4) a similar swelling for large
s. Note finally that the Ginzburg parameter Gz ≈ 1.8 is much larger than for the BFM
systems. As we have emphasized in Sec. II, this should, however, not influence the validity
of the perturbation prediction of the expected 1/
s-swelling of the chains when expressed
in terms of the measured effective bond length.
IV. NUMERICAL RESULTS
As illustrated in Fig. 1, a chain segment of curvilinear length s > 0 is identified by two
monomers n and m = n + s on the same chain. We compute here various moments of
chain segment properties where we ensemble-average over all chains and all start points n.
The statistical accuracy must therefore always decrease for large s. We concentrate first on
the second moment (p = 1) of the segmental size distribution. Higher moments and the
segmental size distribution are discussed in Sec. IVE.
A. The swelling of chain segments
Scale free regime for 1 ≪ s ≪ N . The mean-squared segment size R2e(s) = 〈r2〉 is
presented in the Figs. 4, 5 and 6. The first plot shows clearly that chain segments are swollen,
i.e. R2e(s)/s increases systematically and this up to very large curvilinear distances s. Only
BFM data are shown for clarity. A similar plots exists for the BSM data. In agreement
with Eq. (5) for p = 1, the asymptotic Gaussian behavior (dashed line) is approached
from below and the deviation decays as u(s) ∝ 1/
s (bold line). The bold line indicated
corresponds to be = 3.244 and c1 ≈ ce ≈ 0.41 which fits nicely the data over several decades
in s — provided that chain end effects can be neglected (s ≪ N). Note that a systematic
underestimation of the true effective bond length would be obtained by taking simply the
largest R2e(s)/s ≈ 3.232 value available, say, for monodisperse chains of length N = 2048.
Finite chain-size effects. Interestingly, R2e(s)/s does not approach the asymptotic limit
monotonicly. Especially for short chains one finds a non-monotonic behavior for s → N .
This means that the total chain end-to-end distance Re(s = N − 1) must show even more
pronounced deviations from the asymptotic limit. This is confirmed by the dashed line
representing the b2e(N) ≡ R2e(N − 1)/(N − 1) data points given in Tab. I. We emphasize
that the non-monotonicity of R2e(s)/s becomes weaker with increasing N and that, as one
expects, the inner distances, as well as the total chain size, are characterized by the same
effective bond length be for large s or N . The non-monotonic behavior may be qualitatively
understood by the reduced self-interactions at chain ends which lessens the swelling on these
scales. These finite-N corrections have been calculated analytically using the full Debye
function for the effective interaction potential ṽ(q), Eq. (14). The prediction for the total
chain end-to-end vector given in Eq. (16) is indicated in Fig. 4 (dash-dotted line) where we
have replaced the weakly N -dependent integral I(xu) by its upper bound value for infinite
chains
R2e(N − 1)
b2e (N − 1)
1.59ce√
N − 1
. (19)
We have changed here the chain length N in the analytical formula (obtained for large chains
where N ≈ N − 1) to the curvilinear length N − 1. This is physically reasonable and allows
to take better into account the behavior of small chains. Note that Eq. (19) is similar to
Eq. (5) — apart from a slightly larger prefactor explaining the observed stronger deviations.
Theory compares well with the measured data for large N . It does less so for smaller N ,
as expected, where the chain length dependence of the numerical integral I(xu(N)) ≤ 1.59
must become visible. This explains why the data points are above the dash-dotted line.
Note also that additional non-universal finite-N effects not accounted for by the theory are
likely for small N . In contrast to this, Re(s) is well described by the theory even for rather
small s provided that N is large and chain end effects can be neglected. In summary, it is
clear that one should use the segment size Re(s) rather than the total chain size to obtain
in a computational study a reliable fit of the effective bond length be.
Extrapolation of the effective bond length of asymptotically long chains. The represen-
tation chosen in Fig. 4 is not the most convenient one for an accurate determination of be
and c1. How precise coefficients may be obtained according to Eq. (5) is addressed in the
Figs. 5 and 6. The fitting of the effective bond length be and its accuracy is illustrated in
Fig. 5 for BFM chains of length N = 2048. This may be first done approximately in linear
coordinates by plotting R2e(s)/s as a function of 1/
s (not shown). Since data for large s
are less visible in this representation, we recommend for the fine-tuning of be to switch then
to logarithmic coordinates with a vertical axis y = 1 − R2e(s)/b2es for different trial values
of be. The correct value of be is found by adjusting the vertical axis y such that the data
extrapolates linearly as a function of 1/
s to zero for large s. We assume for the fine-tuning
that higher order perturbation corrections may be neglected, i.e. we take Eq. (5) literally.
(We show below that higher order corrections must indeed be very small.) The plot shows
that this method is very sensitive, yielding a best value that agrees with the theory over
more than one order of magnitude without curvature. As expected, it is not possible to
rationalize the numerically obtained values be ≈ 3.244 for the BFM and be ≈ 1.34 for the
BSM using Eq. (11). According to Eq. (4) these fit values imply the theoretical swelling
coefficients ce = 0.41 for the BFM and ce = 0.44 for the BSM.
Empirical swelling coefficients. As a next step the horizontal axis is rescaled such that
all data sets collapse on the bisection line, i.e. using Eq. (5) we fit for the empirical swelling
coefficient c1 and compare it to the predicted value ce. This rescaling of the axes allows to
compare both models in Fig. 6. For clarity the BSM data have been shifted upwards. For
the BFM we find c1/ce ≈ 1.0, as expected, while our BSM simulations yield a slightly more
pronounced swelling with c1/ce ≈ 1.2.
Segmental radius of gyration. Also indicated in Fig. 6 is the segmental radius of gyration
Rg(s) (filled circles) computed as usual [6] as the variance of the positions of the segment
monomers around their center of mass. Being the sum over all s + 1 monomers, it has a
much better statistics compared to Re(s). The scaling used can be understood by expressing
the radius of gyration R2g(s) =
(s+1)2
l=n r
in terms of displacement vectors rkl
[6]. Using Eq. (5) and integrating twice this yields
6R2g(s)
b2e(s+ 1)
. (20)
Plotting the l.h.s. of this relation against the r.h.s. we obtain a perfect data collapse on the
bisection line where we have used the same parameters be and c1 as for the mean-squared
segment size. This is an important cross-check which we strongly recommend. Different
values indicate insufficient sample equilibration.
B. Chain connectivity and recursion relation
As was emphasized in Sec. IIA the observed swelling is due to an entropic repulsion
between chain segments induced by the incompressibility of the melt. To stress the role of
chain connectivity we repeat the general scaling argument given above in a form originally
proposed by Semenov and Johner for ultrathin films [12]. As shown in Fig. 7 we test the
relation
Kλ(s) ≡
R2e(λs)− λR2e(s)
λ)R2e(s)
≈ u(s) ≡ s
ρRe(s)d
with Kλ(s) being a direct measure of the non-Gaussianity (λ being a positive number)
comparing the size of a segment of length λs with the size of λ segments of length s joined
together. (The prefactor 1/(λ −
λ) in the definition of Kλ(s) has been introduced for
convenience.) Equivalently, this can be read as a measure for the swelling of a chain where
initially the interaction energy u between the segments has been switched off. Kλ is a
functional of u(s) with Kλ[u = 0] = 0. The analytic expansion of the functional must be
dominated by the linear term (as indicated by ≈ in the above relation) simply because u is
very small. Altogether, Eq. (21) yields a recursion relation relating Re(λs) with Re(s) for
any λ provided 1 ≪ s < λs ≪ N . It can be solved, leading (in lowest order) to Eq. (5) with
p = 1. This may be seen from the ansatz R2e(s) = b
es(1− ce/sω−1+ . . .) which readily yields
ω = 3/2 and ce ≈ 1/ρb3e .
Eq. (21) has been validated directly in Fig. 7 for λ = 2 (corresponding to two segments
of length s joined together) for the BFM and the BSM as indicated. In addition, for BFM
chains of length N = 2048 several values of λ have been given. As suggested by Eq. (4),
we have plotted Kλ(s) as a function of (c1/ce) u(s) with u(s) ≡
24/π3s/ρR3e(s) ≈ ce/
The prefactor of u(s) allows a convenient comparison with Fig. 5. Note the perfect data
collapse for all data sets. More importantly, the predicted linearity is well confirmed for
large segments (1 ≪ s) and this without any tunable parameter for the vertical axis, as was
needed in the previous Figs. 5 and 6.
C. Intrachain bond-bond correlations
Expectation from Flory’s hypothesis. An even more striking violation of Flory’s ideality
hypothesis may be obtained by computing the bond-bond correlation function, defined by
the first Legendre polynomial P (s) = 〈lm=n+s · ln〉 /l2 where the average is performed, as
before, over all possible pairs of monomers (n,m = n+ s) [63]. Here, li = ri+1 − ri denotes
the bond vector between two adjacent monomers i and i + 1 and l2 = 〈l2n〉n the mean-
squared bond length. The bond-bond correlation function is generally believed to decrease
exponentially [4]. This belief is based on the few simple single chain models which have
been solved rigorously [4, 40] and on the assumption that all long range interactions are
negligible on distances larger than the screening length ξ. Hence, only correlations along
the backbone of the chains are expected to matter and it is then straightforward to work
out that an exponential cut-off is inevitable due to the multiplicative loss of any information
transferred recursively along the chain [4].
Asymptotic behavior in the melt. That this reasoning must be incorrect follows imme-
diately from the relation
P (s) =
R2e(s) (22)
expressing the bond-bond correlation function as the curvature of the second moment of
the segment size distribution. It is obtained from the identity 〈ln · lm〉 = 〈∂nrn · ∂mrm〉 =
−∂n∂m 〈r2nm〉 /2. (Note that the velocity correlation function is similarly related to the
second derivative of the mean-square displacement with respect to time [41].) Hence, P (s)
allows us to probe directly the non-Gaussian corrections without any ideal contribution.
This relation together with Eq. (5) suggests an algebraical decay P (s) = cP/s
ω with
ω = 3/2 , cP = c1 (be/l)
2/8 ≈
ρl2be
of the bond-bond correlation function for dense solutions and melts, rather than the ex-
ponential cut-off expected from Flory’s hypothesis. This prediction (bold line) is perfectly
confirmed by the larger chains (N > 256) indicated in Fig. 8. In principle, the swelling
coefficient, c1 ∼ cP, may also be obtained from the power law amplitude of the bond-bond
correlation function, however, to lesser accuracy than by the previous method (Fig. 5). One
reason is that P (s) decays very rapidly and does not allow a precise fit beyond s ≈ 102. The
values of cP obtained from c1 are indicated in Tab. II. Data from the BSM have also been
included in the figure to demonstrate the universality of the result. The vertical axis has
been rescaled with cP which allows to collapse the data of both models.
Finite chain-size corrections. As can be seen for N = 16, exponentials are compatible
with the data of short chains. This might explain how the power law scaling has been
overlooked in previous numerical studies, since good statistics for large chains (N > 1000)
has only become available recently. However, it is clearly shown that P (s) approaches
systematically the scale free asymptote with increasing N . The departure from this limit is
fully accounted for by the theory if chain end effects are carefully considered (dashed lines).
Generalizing Eq. (23) and using the Padé approximation, Eq. (17), perturbation theory
yields
P (s) =
1 + 3x+ 5x2
1 + x
(1− x)2 (24)
where we have set x =
s/N . For x ≪ 1 this is consistent with Eq. (23). In the limit of
large s → N , the correlation functions vanish rigorously as P (s) ∝ (1 − x)2. Considering
that non-universal features cannot be neglected for short chain properties and that the
theory does not allow for any free fitting parameter, the agreement found in Fig. 8 is rather
satisfactory.
D. Higher moments and associated coefficients
Effective bond length and empirical swelling coefficients. The preceding discussion fo-
cused on the second moment of the segmental size distribution G(r, s). We have also com-
puted for both models higher moments 〈r2p〉 with p ≤ 5. If traced in log-linear coordinates
as y = (6pp! 〈r2p〉 /(2p + 1)!sp)1/p vs. x = s higher moments approach b2e from below —
just as the second moment presented in Fig. 4. The deviations from ideality are now more
pronounced and increase with p (not shown). The moments are compared in Fig. 6 with
Eq. (5) where they are rescaled as y = Kp(s) as defined in Eq. (1) and plotted as functions
of x =
3(2pp!p)2
2(2p+1)!
. The prediction is indicated by bold lines. It is important that the same
effective bond length be is obtained from the analysis of all functions Kp(s) as illustrated in
Fig. 6. Otherwise we would regard equilibration and statistics as insufficient.
The empirical swelling coefficients cp are obtained, as above in Sec. IVA, by shifting the
data horizontally. A good agreement with the expected cp/ce ≈ 1 is found for both models
and all moments as may be seen from Tab. II. This confirms the renormalization of the
Kuhn segment b → be of the Gaussian reference chain in agreement with our discussion in
Sec. II B. Otherwise we would have measured empirical coefficients decreasing strongly as
cp/ce ≈ (b/be)2p−3 with p. Since the effective bond length of non-interacting chains are known
for the BFM (b ≈ 2.688) and the BSM (b ≈ 0.97), one can simply check, say for p = 5, that
the non-renormalized values would correspond to the ratios c5/ce ≈ (2.688/3.244)7 ≈ 0.3 for
the BFM and c5/ce ≈ (0.97/1.34)7 ≈ 0.1 for the BSM. This is clearly not consistent with
our data.
It should be emphasized that both coefficients be and cp are more difficult to determine
for large p, since the linear regime for x ≪ 1 in the representation chosen in Fig. 6 becomes
reduced. For large x ≫ 1 one finds that y(x) → 1, i.e. 〈r2p〉 /b2pe sp → 0. This trivial
departure from both Gaussianity and the 1/
s-deviations we try to describe, is due to the
finite extensibility of chain segments of length s which becomes more marked for larger
moments probing larger segment sizes. The data collapse for both x-regimes is remarkable,
however. Incidentally, it should be noted that for the BSM the empirical swelling coefficients
are slightly larger than expected. At present we do not have a satisfactory explanation for
this altogether minor effect, but it might be attributed to the fact that neighbouring BSM
beads along the chain strongly interpenetrate — an effect not considered by the theory.
Non-Gaussian parameter αp. The failure of Flory’s hypothesis can also be demonstrated
by means of the standard non-Gaussian parameter
αp(s) ≡ 1−
(2p+ 1)!
〈r2pnm〉
〈r2nm〉
p (25)
comparing the 2p-th moment with the second moment (p = 1). In contrast to the closely
related parameter Kp(s) this has the advantage that here two measured properties are com-
pared without any tuneable parameter, such as be, which has to be fitted first. Fig. 9 presents
αp(s) vs. ce/
s for the three moments with p = 2, 3, 4. For each p we find perfect data col-
lapse for all chain lengths and both models and confirm the linear relationship αp(s) ≈ u(s)
expected. The lines indicate the theoretical prediction
αp(s) =
3 (2pp!p)2
2 (2p+ 1)!
which can be derived from Eq. (5) by expanding the second moment in the denominator. An
alternative derivation based on the coefficients of the expansion of the generating function
G(q, s) in q2 is indicated by Eq. (A2) in the Appendix. Having confirmed above that cp/ce ≈
1, we assume in Eq. (26) that cp = ce to simplify the notation. The prefactors 6/5, 111/35
and 604/105 for p = 2, 3 and 4 respectively are nicely confirmed. They increase strongly
with p, i.e. the non-Gaussianity becomes more pronounced for larger moments as already
mentioned. Note also the curvature of the data at small s due to the finite extensibility
of the segments which becomes more marked for higher moments. If one plots αp(s) as a
function of the r.h.s. of Eq. (26) all data points for all moments and even for too small s
collapse on one master curve (not shown) — just as we have seen before in Fig. (6).
Correlations of different directions. A similar correlation function is presented in Fig. 10
which measures the non-Gaussian correlations of different spatial directions. It is defined by
Kxy(s) ≡ 1−
〈x2 y2〉
〈x2〉 〈y2〉
for the two spatial components x and y of the vector r as illustrated by the sketch given
at the bottom of Fig. 10. Symmetry allows to average over the three pairs of directions
(x, y), (x, z) and (x, z). Following the general scaling argument given in Sec. II we expect
Kxy(s) ≈ u(s) ≈ ce/
s which is confirmed by the perturbation result
Kxy(s) =
= K2(s). (28)
This is nicely confirmed by the linear relationship found (bold line) on which all data from
both simulation models collapse perfectly. The different directions of chain segments are
therefore coupled. As explained in the Appendix (Eq. (A3)), Kxy(s) and α2(s) must be
identical if the Fourier transformed segmental size distribution G(q, s) can be expanded
in terms of q2 and this irrespective of the values the expansion coefficients take. Fig. 10
confirms, hence, that our computational systems are perfectly isotropic and tests the validity
of the general analytical expansion.
The correlation function Kxy is of particular interest since the zero-shear viscosity should
be proportional to
∼ 〈x2y2〉 = 〈x2〉 〈y2〉 (1 − Kxy(s)). We assume here following
Edwards [6] that only intrachain stresses contribute to the shear stress σxy. Hence, our
results suggest that the classical calculations [6] — assuming incorrectly Kxy = 0 — should
be revisited.
E. The segmental size distribution
We turn finally to the segmental size distribution G(r, s) itself which is presented in
Figs. 11, 12 and 13. From the theoretical point of view G(r, s) is the most fundamental
property from which all others can be derived. It is presented last since it is computationally
more demanding — at least if high accuracy is needed — and coefficients such as be may
be best determined directly from the moments. The normalized histograms G(r, s) are
computed by counting the number of segment vectors between r − dr/2 and r + dr/2 with
dr being the width of the bin and one divides then by the spherical bin volume. Since
the BFM model is a lattice model, this volume is not 4πr2dr but given by the number of
lattice sites the segment vector can actually point to for being allocated to the bin. Incorrect
histograms are obtained for small r if this is not taken into account. (Averages are taken
over all segments and chains, just as before.) Clearly, non-universal physics must show up
for small vector length r and small curvilinear distance s and we concentrate therefore on
values r ≫ σ and s ≥ 31.
When plotted in linear coordinates as in Fig. 11, G(r, s) compares roughly with the
Gaussian prediction G0(r, s) given by Eq. (2), but presents a distinct depletion for small
segment sizes with n ≡ r/be
s ≪ 1 and an enhanced regime for n ≈ 1. A second depletion
region for large n ≫ 1 — expected from the finite extensibility of the segments — can be best
seen in the log-log representation of the data (not shown). To analyse the data it is better
to consider instead of G(r, s) the relative deviation δG(r, s)/G0(r, s) = G(r, s)/G0(r, s)− 1
which should further be divided by the strength of the segmental correlation hole, ce/
As presented in Fig. 12 this yields a direct test of the key relation Eq. (3) announced in the
Introduction. The figure demonstrates nicely the scaling of the data for all s and for both
models. It shows further a good collapse of the data close to the universal function f(n)
predicted by theory (bold line). Note that the depletion scales as 1/n for small segment
sizes (dashed line). The agreement of simulation and theory is by all standards remarkable.
(Obviously, error bars increase strongly for n ≫ 1 where G0(r, s) decreases strongly. The
regime for very large n where the finite extensibility of segments matters has been omitted
for clarity.) We emphasize that this scaling plot depends very strongly on the value be which
is used to calculate the Gaussian reference distribution.
If a precise value is not available we recommend to use instead the scaling variable
m = r/Re(s) for the horizontal axis, i.e. to replace the scale be
s estimated from the
behaviour of asymptotically long chains by the measured (mean-squared) segment size
for the given s. The Gaussian reference distribution is then accordingly G0(m,Re(s)) =
(3/2πRe(s)
2) exp(−3
m2). The corresponding scaling plot is given in Fig. 13. It is simi-
lar and of comparable quality as the previous plot. Changing the scaling variable from
n = r/be
s to m = r/Re(s) ≈ (r/be
s)(1 + ce/2
s) changes somewhat the universal func-
tion. Expanding the previous result, Eq. (3), this adds even powers of m to the function
f(n) given in Eq. (3)
f(n) ⇒ f(m) =
+ 9m+
m2 − 9
. (29)
That the two additional terms in the function are correct can be seen by computing the
second moment 4π
drr4δG(r, s) which must vanish by construction. The rescaled relative
deviation is somewhat broader than in the previous plot due to the additional term scaling
as m2. As already stressed this scaling does not rely on the effective bond length be and
is therefore more robust. It has the nice feature that it underlines that there is only one
characteristic length scale relevant for the swelling induced by the segmental correlation
hole, the typical size of the chain segment itself.
V. CONCLUSION
Issues covered and central theoretical claims. We have revisited Flory’s famous ideality
hypothesis for long polymers in the melt by analyzing both analytically and numerically
the segmental size distribution G(r, s) and its moments for chain segments of curvilinear
length s. We have first identified the general mechanism that gives rise to deviations from
ideal chain behavior in dense polymer solutions and melts (Sec. II). This mechanism rests
upon the interplay of chain connectivity and the incompressibility of the system which
generates an effective repulsion between chain segments (Fig. 2). This repulsion scales like
u(s) ≈ ce/
s where the “swelling coefficient” ce ≈ 1/b3eρ sets the strength of the interaction.
It is strong for small segment length s, but becomes weak for s → N in the large-N limit.
The overall size of a long chain thus remains almost ‘ideal’, whereas subchains are swollen as
described by Eq. (5). Most notably, the relative deviation δG(r, s)/G0(r, s) of the segmental
size distribution from Gaussianity should be proportional to u(s). As a function of segment
size r, the repulsion manifests itself by a strong 1/r-depletion at short distances r ≪ be
and a subsequent shift of the histogram to larger distances (Eq. (3)).
Summary of computational results. Using Monte Carlo and molecular dynamics simu-
lation of two coarse-grained polymer models we have verified numerically the theoretical
predictions for long and flexible polymers in the bulk. We have explicitly checked (e.g.,
Figs. 7, 9, 13) that the relative deviations from Flory’s hypothesis scale indeed as 1/
Especially, the measurement of the bond-bond correlation function P (s), being the second
derivative of the second moment of G(r, s) with respect of s, allows a very precise verifi-
cation (Fig. 8) and shows that higher order corrections beyond the first-order perturbation
approximation must be small. The most central and highly non-trivial numerical verification
concerns the data collapse presented in Figs. 12 and 13 for the segmental size distribution
of both computational models. All other statements made in this paper can be derived and
understood from this key finding. It shows especially that the swelling coefficient ce must
be close to the predicted value, Eq. (4).
It is well known [10] that the effective bond length is difficult to predict at low com-
pressibility and no attempt has been done to do so in this paper. We show instead how the
systematic swelling of chain segments – once understood – may be used to extrapolate for
the effective bond length of asymptotically long chains. Figs. 5 and 6 indicate how this may
be done using Eq. (5). The high precision of our data is demonstrated in Fig. 12 by the
successful scaling of the segmental size distribution.
For several moments 〈r2p〉 we have also fitted empirical swelling coefficients cp using
Eq. (5). In contrast to the effective bond length be these coefficients are rather well pre-
dicted by one-loop perturbation theory if the bond length b of the reference Hamiltonian is
renormalized to the effective bond length be, as we have conjectured in Sec. II B. Since the
empirical swelling coefficients, cp ≈ ce(b/be)2p−3, would otherwise strongly depend on the
moment taken, as shown in Eq. (12), our numerical data (Tab. II) clearly imply b/be ≈ 1.
Minor deviations found for the BSM samples may be attributed to the fact that monomers
along the BSM chains do strongly overlap — an effect not taken into account by the theory.
To clarify ultimately this issue we are currently performing a numerical study where we
systematically vary both the compressibility and the bond length of the BSM.
General background and outlook. The most striking result presented in this work con-
cerns the power law decay found for the bond-bond correlation function, P (s) ∝ 1/s3/2
(Fig. 8). This result suggests an analogy with the well-known long-range velocity correla-
tions found in dense fluids by Alder and Wainwright nearly fourty years ago [41, 42]. In
both cases, the ideal uncorrelated object is a random walker which is weakly perturbed (for
d > 2) by the self-interactions generated by global constraints. Although these constraints
are different (momentum conservation for the fluid, incompressibility for polymer melts) the
weight with which these constraints increase the stiffness of the random walker is always
proportional to the return probability. It can be shown that the correspondence of both
problems is mathematically rigorous if the fluid dynamics is described on the level of the
linearized Navier-Stokes equations [43].
We point out that the physical mechanism which has been sketched above is rather
general and should not be altered by details such as a finite persistence length — at least
not as long as nematic ordering remains negligible and the polymer chains are sufficiently
long. (Similarly, velocity correlations in dense liquids must show an analytical decay for
sufficiently large times irrespective of the particle mass and the local static structure of
the solution.) While this paper focused exclusively on scales beyond the correlation length
of the density fluctuations, i.e. qξ ≪ 1 or s/g ≫ 1, where the polymer solution appears
incompressible, effects of finite density and compressibility can be readily described within
the same theoretical framework and will be presented elsewhere [43]. To test our predictions,
flexible chains should be studied preferentially, since the chain length required for a clear-
cut description increases strongly with persistence length. This is in fact confirmed by
preliminary and on-going simulations using the BSM algorithm.
In this work we have only discussed properties in real space as a function of the curvilin-
ear distance s. These quantities are straightforward to compute in a computer simulation
but are barely experimentally relevant. The non-Gaussian deviations induced by the seg-
mental correlation hole arise, however, also for an experimentally accessible property, the
intramolecular form factor (single chain scattering function) F (q). As explained at the
end of the Appendix, the form factor can be readily obtained by integrating the Fourier
transformed segmental size distribution given in Eq. (3). This yields
q2F (q) ≈ 12
in agreement with the result obtained in Refs. [14, 19] by direct calculation of the form factor
for very long equilibrium polymers. As a consequence of this, the Kratky plot (q2F (q) vs.
wave vector q) should not exhibit the plateau expected for Gaussian chains in the scale-free
regime, but rather noticeable non-monotonic deviations. See Fig. 3 of [19]. This result
suggests to revisit experimentally this old pivotal problem of polymer science.
Our work is part of a broader attempt to describe systematically the effects of correlated
density fluctuations in dense polymer systems, both for static [12, 13, 44, 45] and dynam-
ical [29, 35, 46] properties. An important unresolved question is for instance whether the
predicted long-range repulsive forces of van der Waals type (“Anti-Casimir effect”) [13, 45]
are observable, for instance in the oscillatory decay of the standard density pair-correlation
function of dense polymer solutions. Since the results presented here challenge an important
concept of polymer physics, they should hopefully be useful for a broad range of theoreti-
cal approaches which commonly assume the validity of the Gaussian chain model down to
molecular scales [47, 48, 49]. This study shows that a polymer in dense solutions should not
be viewed as one soft sphere (or ellipsoid) [50, 51, 52], but as a hierarchy of nested segmental
correlation holes of all sizes aligned and correlated along the chain backbone (Fig. 2 (b)). We
note that similar deviations from Flory’s hypothesis have been reported recently for linear
polymers [16, 17, 47] and polymer gels and networks [53, 54]. The repulsive interactions
should also influence the polymer dynamics, since strong deviations from Gaussianity are
expected on the scale where entanglements become important, hence, quantitative predic-
tions for the entanglement length Ne have to be regarded with more care. The demonstrated
swelling of chains should be included in the popular primitive path analysis for obtaining Ne
[55], especially if ‘short’ chains (N < 500) are considered. The effect could be responsible
for observed deviations from Rouse behavior [26, 56] as may be seen by considering the cor-
relation function Cpq ≡ 〈Xp ·Xq〉 of the Rouse modes Xp = 1N
dnrn cos(npπ/N) where
p, q = 0, . . . , N − 1 [6, 57]. Using (rn − rm)2 = r2n + r2m − 2rn · rm for the segment size, this
correlation function can be readily expressed as an integral over the second moment of the
segmental size distribution
Cpq = −
(rn − rm)2
cos(npπ/N) cos(mpπ/N) (31)
which can be solved using our result Eq. (5). This implies for instance for p = q that
Cpp =
2(πp)2
1− π√
. (32)
The bracket entails an important correction with respect to the classical description given by
the prefactor [6]. We are currently working out how static corrections, such as those for Cpp,
may influence the dynamics for polymer chains without topological constraints. (This may
be realized, e.g., within the BFM algorithm by using the L26 moves described in Sec. IIIA.)
Moreover, for thin polymer films of width H the repulsive interactions are known to be
stronger than in the bulk [12]. This provides a mechanism to rationalize the trend towards
swelling observed experimentally [58] and confirmed computationally [21]:
= log(s)/H. (33)
(Prefactors omitted for clarity.) Here Rx(s) and bx denote the components of the segment
size and the effective bond length parallel to the film. It also explains the (at first sight
surprising) systematic increase of the polymer dynamics with decreasing film thickness [22].
Specifically, the parallel component of the monomer mean-squared displacement gx(t) is
expected to scale as gx(t) ≈ R2x(s(t)) ∝ t1/4(1 + log(t)/H) for long reptating chains where
s(t) ∝ t1/4 [6]. (The corresponding effect for the three-dimensional bulk should be small,
however.) For the same reason (flexible) polymer chains close to container walls must be
more swollen and, hence, faster on intermediate time scales than their peers in the bulk.
Acknowledgments
We thank T. Kreer, S. Peter and A.N. Semenov (all ICS, Strasbourg, France),
S.P. Obukhov (Gainesville, Florida) and M. Müller (Göttingen, Germany) for helpful discus-
sions. A generous grant of computer time by the IDRIS (Orsay) is also gratefully acknowl-
edged. J.B. acknowledges financial support by the IUF and from the European Community’s
“Marie-Curie Actions” under contract MRTN-CT-2004-504052.
APPENDIX A: MOMENTS OF THE SEGMENTAL SIZE DISTRIBUTION AND
THEIR GENERATING FUNCTION
Higher moments of the segmental size distribution G(r, s) can be systematically obtained
from its Fourier transformation
G(q, s) =
d3r G(r, s) exp(iq · r),
which is in this context sometimes called the “generating function” [59]. For an ideal
Gaussian chain, the generating function is then G0(q, s) = exp(−sq2a2) where we have
used a2 = b2/6 instead of the bond length b2 to simplify the notation. Moments of the
size distribution are given by proper derivatives of G(q, s) taken at q = 0. For example,
〈r2p〉 = (−1)p∆pG(q, s)q=0 (with ∆ being the Laplace operator with respect to the wave vec-
tor q). A moment of order 2p is, hence, linked to only one coefficient A2p in the systematic
expansion, G(q, s) =
p=0A2pq
2p, of G(q, s) around q = 0. For our example this implies
= (−1)p(2p+ 1)! A2p (A1)
in general and more specifically for a Gaussian distribution 〈r2p〉0 =
(2p+1)!
spa2p. The non-
Gaussian parameters read, hence,
αp(s) ≡ 1−
(2p+ 1)!
〈r2p〉
〈r2〉p
= 1− p! A2p
, (A2)
which implies (by construction) αp = 0 for a Gaussian distribution. As various moments of
the same global order 2p are linked to the same A2p they differ by a multiplicative constant
independent of the details of the (isotropic) distribution G(q, s). For example, 〈r2〉 = 6|A2|,
〈r4〉 = 120A4, 〈x2〉 = 〈y2〉 = 2|A2|, 〈x2y2〉 = 8A4 with x and y denoting the spatial
components of the segment vector r. Using Eq. (A2) for p = 2 it follows that
Kxy(s) ≡ 1−
〈x2y2〉
〈x2〉〈y2〉
= 1− 2A4
= α2(s), (A3)
i.e. the properties α2(s) and Kxy(s) discussed in Figs. 9 and 10 must be identical in general
provided that G(q, s) is isotropic and can be expanded in q2.
We turn now to specific properties ofG(q, s) computed for formally infinite polymer chains
in the melt. In practice, these results are also relevant for small segments in large chains,
N ≫ s ≫ 1, and, especially, for segments located far from the chain ends. These chains are
nearly Gaussian and the generating function can be written as G(q, s) = G0(q, s) + δG(q, s)
where δG(q, s) = −〈UG〉0 + 〈U〉0〈G〉0 is a small perturbation under the effective interaction
potential ṽ(q) given by Eq. (9). To compute the different integrals it is more convenient to
work in Fourier-Laplace space (q, t) with t being the Laplace variable conjugate to s:
δG(q, t) =
ds δG(q, s)e−st.
As illustrated in Fig. 14, there a three contributions to this perturbation: one due to in-
teractions between two monomers inside the segment (left panel), one due to interactions
between an internal monomer and an external one (middle panel) and one due to interac-
tions between two external monomers located on opposite sides (right panel). In analogy to
the derivation of the form factor described in Ref. [14] this yields:
δG(q, t) = −
(q2a2 + t)2
q2a2 + t−
(q2a2 + t)2
4πqa2ξ2
Arctan
a/ξ +
a/ξ +
q2a2 + t
q2a2 + t
4πqa4
Arctan
a/ξ +
a/ξ +
q2a2 + t
4πqa6
Arctan
q2a2 + t
. (A4)
The graph given in the left panel of Fig. 14 corresponds to the first two lines, the middle
panel to the third line and the right panel to the last one. Seeking for the moments we
expand δG(q, t) around q = 0. Having in mind chain strands counting many monomers
(s ≫ 1), we need only to retain the most singular terms for t → 0. Defining the two
dimensionless constants d = vξ/3πa4 = 12vξ/πb4 and c = (3π3/2a3ρ)−1 =
24/π3/b3ρ this
expansion can be written as
δG(q, t) = −
d a2q2 +
Γ(3/2)
c a2q2 + . . . (A5)
d a4q4 − 1
Γ(5/2)
c a4q4 + . . .
d a6q6 +
Γ(7/2)
c a6q6 + . . .
+ . . .
where we have used Euler’s Gamma function Γ(α) [60]. The first leading term at each order
in q2 — being proportional to the coefficient d — ensures the renormalization of the effective
bond length. The next term scaling with the coefficient c corresponds to the leading finite
strand size correction. Performing the inverse Laplace transformation Γ(α)/tα → sα−1 and
adding the Gaussian reference distribution G0(q, s) this yields the A2p-coefficients for the
expansion of G(q, s) around q = 0:
A0 = 1
A2 = −a2s
1 + d− c√
1 + 2d−
A6 = −
1 + 3d− 216
A8 = . . . (A6)
More generally, one finds
A2p =
(−1)p
(sa2)p
1 + pd− 3(2
pp!p)2
2(2p+ 1)!
From this result and using Eq. (A1) one immediately verifies that the moments of the
distribution are given by the Eqs. (11) and (12). Using Eq. (A2) one justifies similarly
Eq. (26) for the non-Gaussian parameter αp.
These moments completely determine the segmental distribution G(r, s) which is indi-
cated in Eq. (13). While at least in principle this may be done directly by inverse Fourier-
Laplace transformation of the correction δG(q, t) to the generating function it is helpful to
simplify further Eq. (A4). We observe first that δG(q, t) does diverge for strictly incom-
pressible systems (v → ∞) and one must keep v finite in the effective potential whenever
necessary to ensure convergence (actually everywhere but in the diagram corresponding to
the interaction between two external monomers). Since we are not interested in the wave
vectors larger than 1/ξ we expand δG(q, t) for ξ → 0 which leads to the much simpler
expression
δG(q, t) ≈ − vξq
3πa2(a2q2 + t)2
t(3a2q2 + t)
(a2q2 + t)2
Arctan[ aq√
+ o(vξ3). (A8)
The first term diverges as
v for diverging v. It renormalizes the effective bond length in
the zero order term which is indicated in the first line of Eq. (13). The next two terms scale
both as v0. Subsequent terms must all vanish for diverging v and can be discarded. It is
then easy to perform an inverse Fourier-Laplace transformation of the two relevant v0 terms.
This yields
δG(x, s) = G0(x, s)
with x = r/a
6n. This is consistent with the expression given in the second line of
Eq. (13).
We note finally that the intramolecular form factor F (q) = 1
n,m=1 〈exp(iq · (rn − rm)〉
of asymptotically long chains can be readily obtained from Eq. (A8). Observing that
〈exp(iq · (rn − rm)〉 =
d3r exp(iq · r)G(r, s) = G(q, s) one finds
δF (q) = 2
ds δG(q, s) = 2 δG(q, t = 0) = −2 vξ
, (A10)
where we used the third term of Eq. (A8) in the last step. The first term in Eq. (A8) is
discarded as before, since it renormalizes the effective bond length in the reference form
factor: F0(q) = 12/b
2q2 ⇒ 12/b2eq2. It follows, hence, that within first-order perturbation
theory
F (q) = F0(q) + δF (q) ≈ F0(q)
(A11)
as indicated by Eq. (30) in the Conclusion. This is equivalent to the result 1/F (q)−1/F0(q) ≈
q3/32ρ discussed in Refs. [14, 19] for polymer melts and anticipated by Schäfer [11] by
renormalization group calculations of semidilute solutions.
N nch τe Re Rg be(N)
6bg(N)
16 216 1214 11.7 4.8 2.998 2.939
32 215 3485 17.1 7.0 3.066 3.030
64 214 1.1 · 104 24.8 10.1 3.116 3.094
128 8192 3.3 · 104 35.6 14.5 3.153 3.139
256 4096 1.0 · 105 50.8 20.7 3.179 3.171
512 2048 3.2 · 105 72.2 29.5 3.200 3.193
1024 1024 1.0 · 106 103 42.0 3.216 3.212
2048 512 3.2 · 106 146 59.5 3.227 3.223
4096 256 9.7 · 106 207 85.0 3.235 3.253
8192 128 2.9 · 107 294 120 3.249 3.248
TABLE I: Various static properties of dense BFM melts of number density ρ = 0.5/8: the chain
length N , the number of chains nch per box, the relaxation time τe characterized by the diffusion
of the monomers over the end-to-end distance and corresponding to the circles indicated in Fig. 3,
the root-mean-squared chain end-to-end distance Re and the radius of gyration Rg of the total
chain (s = N − 1). The last two columns give estimates for the effective bond length from the
end-to-end distance, be(N) ≡ Re/(N − 1)1/2, and the radius of gyration, bg(N) ≡ Rg/
N . The
dashed line in Fig. 4 indicates be(N)
2. Apparently, both estimates increase monotonicly with N
reaching be(N) ≈
6bg(N) ≈ 3.2 for the largest chains available. Note that
6bg(N) < be(N) for
smaller N .
Property BFM BSM
Length unit lattice constant bead diameter
Temperature kBT 1 1
Number density ρ 0.5/8 0.84
Linear box size L 256 ≤ 62
Number of monomers nmon 1048576 ≤ 196608
Largest chain length N 8192 1024
Mean bond length 〈|ln|〉 2.604 0.97
l = 〈l2n〉1/2 2.636 0.97
Effective bond length be 3.244 1.34
ρb3e 2.13 2.02
C∞ = (be/l)
2 1.52 1.91
lp = l(C∞ + 1)/2 3.32 1.41
24/π3/ρb3e 0.41 0.44
c1/ce 1.0 1.2
c2/ce 1.0 1.1
c3/ce 1.0 1.0
c4/ce 1.1 1.2
c5/ce 1.1 0.9
cP = c1(be/l)
2/8 0.078 0.124
Dimensionless compressibility g 0.245 0.08
Compression modulus v ≡ 1/gρ 66.7 14.9
vρ/b3eρ 0.96 1.8
TABLE II: Comparison of some static properties of dense BFM and BSM melts. The first six rows
indicate conventions and operational parameters. The effective bond length be and the swelling
coefficients cp (defined in Eq. (5)) are determined from the first five even moments of the segmental
size distribution. The dimensionless compressibility g = S(q → 0)/ρ has been obtained from the
total static structure factor S(q) = 1
∑nmon
k,l=1 〈exp(iq · (rk − rl))〉 in the zero wave vector limit
as shown at the end of Ref. [14]. The values indicated correspond to the asymptotic long chain
behavior. Properties of very small chains deviate slightly.
I ~ −3 I ~ −9
I ~ 45
I ~ −9−
r   n m
FIG. 1: (Color online) Sketch of a polymer chain of length N in a dense melt in d = 3 dimensions.
As notations we use ri for the position vector of a monomer i, li = ri+1 − ri for its bond vector,
r = rm−rn for the end-to-end vector of the chain segment between the monomers n and m = n+s
and r = ||r|| for its length. Segment properties, such as the 2p-th moments
, are averaged over
all possible pairs of monomers (n,m) of a chain and over all chains. The second moment (p = 1)
is denoted Re(s) =
, the total chain end-to-end distance is Re(s = N − 1). The dashed
lines show the relevant graphs of the analytical perturbation calculation outlined in Sec. II B.
The numerical factors indicate for infinite chains (without chain end effects) the relative weights
contributing to the 1/
s-swelling of Re(s) indicated in Eq. (10).
a bdensity ρ
c(r,s)
correlation hole
Segmental
= const
R(s) Repulsion
FIG. 2: (Color online) Role of incompressibility and chain connectivity in dense polymer solutions
and melts. (a) Sketch of the segmental correlation hole of a marked chain segment of curvilin-
ear length s. Density fluctuations of chain segments must be correlated, since the total density
fluctuations (dashed line) are small. Consequently, a second chain segment feels an entropic re-
pulsion when both correlation holes start to overlap. (b) Self-similar pattern of nested segmental
correlation holes of decreasing strength u(s) ≈ s/ρR(s)3 ≈ ce/
s aligned along the backbone of
a reference chain. The large dashed circle represents the classical correlation hole of the total
chain (s ≈ N) [5]. This is the input of recent approaches to model polymer chains as soft spheres
[50, 52]. We argue that incompressibility on all scales and chain connectivity leads to a short
distance repulsion of the segmental correlation holes, which increases with decreasing s.
L06 (conserved topology)
L26,SS
L26,SS,DB
FIG. 3: (Color online) Diffusion time τe over the (root-mean-squared) chain end-to-end distance
Re(N − 1) as a function of chain length N for different versions of the Bond Fluctuation Model
(BFM). All data indicated are for the high number density (ρ = 0.5/8) corresponding to a polymer
melt with half the lattice sites being occupied. We have obtained τe = R
e(N−1)/6Ds from the self-
diffusion coefficient Ds measured from the free diffusion limit of the mean-squared displacement
of all monomers
δr(t)2
= 6Dst. Data from the classical BFM with topology conserving local
Monte Carlo (MC) moves in 6 spatial directions (L06) [26] are represented as stars. All other data
sets use topology violating local MC moves in 26 lattice directions (L26). If only local moves are
used, L26-dynamics is even at relatively short times perfectly Rouse like which allows the accurate
determination of Ds although the monomers possibly have not yet moved over Re(N − 1) for the
largest chain lengths considered. Additional slithering snake (SS) moves increase the efficiency of
the algorithm by approximately an order of magnitude (squares,bold line). The power law exponent
is changed from 2 to an empirical 1.62 (dashed line) if in addition we perform double-bridging (DB)
moves.
N=128
N=256
N=512
N=1024
N=2048
N=4096
N=8192
=3.244
Eq.(5)
FIG. 4: (Color online) Mean-squared segment size Re(s)
2/s vs. curvilinear distance s. We present
BFM data for different chain length N at number density ρ = 0.5/8. The averages are taken over
all possible monomer pairs (n,m = n+s). The statistics deteriorates, hence, for large s. Log-linear
coordinates are used to emphasize the power law swelling over several orders of magnitude of s.
The data approach the asymptotic limit (horizontal line) from below, i.e. the chains are swollen.
This behavior is well fitted by Eq. (5) for 1 ≪ s ≪ N (bold line). Non-monotonic behavior is found
for s → N , especially for small N . The dashed line indicates the measured total chain end-to-end
distances, be(N)
2 ≡ Re(N − 1)2/(N − 1) from Tab. I, showing even more pronounced deviations
from the asymptotic limit. The dash-dotted line compares this data with Eq. (19).
=3.235
=3.240
=3.244
=3.250
=3.255
BFM ρ=0.5/8, N=2048: 
=3.244, c
=0.41=c
FIG. 5: (Color online) Replot of the mean-squared segment size as y = K1(s) = 1 − Re(s)2/b2es
vs. x = c1/
s, as suggested by Eq. (5), for different trial effective bond lengths be as indicated.
Only BFM chains of length N = 2048 are considered for clarity. This procedure is very sensitive
to the value chosen and allows for a precise determination. It assumes, however, that higher order
terms in the expansion of K1(s) may be neglected. The value be is confirmed from a similar test
for higher moments (Fig. 6).
x ~1/s
BFM b
=3.244
BSM b
=1.34
too small s !
FIG. 6: (Color online) Critial test of Eq. (5) where the rescaled moments y = Kp(s) of the segment
size distribution (defined in Eq. (1)) are plotted vs. x =
3(2pp!p)2
2(2p+1)!
. We consider the first five even
moments (p = 1, . . . , 5) for the BFM with N = 2048 and the BSM with N = 1024. Also indicated
is the rescaled radius of gyration, y = 5/8 (1 − 6R2g(s)/b2e(s + 1)), as a function of x = c1/
(filled circles). The BSM data has been shifted upwards for clarity. Without this shift a perfect
data collapse is found for both models and all moments. Keeping the same effective bond length
be for all moments of each model we fit for the swelling coefficients cp by rescaling the horizontal
axis. We find be ≈ 3.244 for the BFM and 1.34 for the BSM. If be is chosen correctly, all data
sets extrapolate linearly to zero for large s (x → 0). The swelling coefficients found are close the
theoretical prediction ce, as indicated in Tab. II. The plot demonstrates that the non-Gaussian
deviations scale as the segmental correlation hole, u(s) ∼ ce/
s and this for all moments as long
as x ≪ 1. The saturation at large x is due to the finite extensibility of short chain segments. Since
this effect becomes more marked for larger moments, the fit of be is best performed for p = 1.
u(s)  c
Kλ(s)
λ=2: N=1024
λ=2: N=2048
λ=2: N=4096
λ=4: N=2048
λ=8: N=2048
λ=16: N=2048
Insufficient statistics & Finite chain length effects
⇒ ω = 3/2
FIG. 7: (Color online) Plot of Kλ(s) as a function of u(s)c1/ce ∼ 1/
s using the measured
u(s) ≡
24/π3s/ρRe(s)
3. For λ = 2 (corresponding to two segments being connected) BFM and
BSM data are compared. Several λ values are given for N = 2048 BFM chains. For chain segments
with 1 ≪ s ≪ N all data sets collapse on the bisection line confirming the so-called “recursion
relation” Kλ ≈ u proposed by Semenov and Johner [12]. The statistics becomes insufficient for
large s (left bottom corner). Systematic deviations arise for s → N due to additional finite-N
effects.
N=128
N=512
N=2048
N=4096
N=8192
Slope ω=3/2
Chain end effects
FIG. 8: (Color online) The bond-bond correlation function P (s)/cP as a function of the curvilinear
distance s. Various chain lengths are given for BFM. Provided that 1 ≪ s ≪ N , all data sets
collapse on the power law slope with exponent ω = 3/2 (bold line) as predicted by Eq. (23). The
dash-dotted curve P (s) ≈ exp(−s/1.5) shows that exponential behavior is only compatible with
very small chain lengths. The dashed lines correspond to the theoretical prediction, Eq. (24), for
short chains with N = 16, 32, 64 and 128 (from left to right).
BFM N=1024
BFM N=2048
BFM N=4096
BFM N=8192
BSM N=1024
y = 6x/5
y = 111x/35
y = 605x/105
small s !
Noise !
FIG. 9: (Color online) Non-Gaussian parameter αp(s) computed for the end-to-end distance of
chain segments as a function of ce/
s. Perfect data collapse for all chain lengths and both sim-
ulation models is obtained for each p. A linear relationship over nearly two orders of magnitude
is found as theoretically expected. Data for three moments (p = 2, 3, 4) are indicated showing
a systematic increase of non-Gaussianity with p. The data curvature for small s becomes more
pronounced for larger p.
N=256
N=512
N=1024
N=2048
N=4096
N=8192
m=n+s
FIG. 10: (Color online) Plot of Kxy(s) = 1−
averaged over all pairs of monomers
(n,m = n+s) and three different direction pairs as a function of ce/
s. As indicated by the sketch
at the bottom of the figure, Kxy(s) measures the correlation of the components of the segment
vector r. All data points collapse and show again a linear relationship Kxy ≈ u(s). Different
directions are therefore coupled! No curvature is observed over two orders of magnitude confirming
that higher order perturbation corrections are negligible. Noise cannot be neglected for large
s > 100 and finite segment-size effects are visible for s ≈ 1.
0.0 0.5 1.0 1.5 2.0
n=r/b
BFM s=32
BFM s=64
BFM s=128
BSM s=32
BSM s=64
BSM s=128
GaussEnhancem
entDepletion
FIG. 11: (Color online) Segment size distribution y = G(r, s)(bes
1/2)3 vs. n = r/bes
1/2 for several
s as indicated in the figure. Only data for BFM with N = 2048 and BSM with N = 1024 are
presented. (A similar plot can be achieved by renormalizing the axes using Re(s) instead of bes
1/2).
The bold line denotes the Gaussian behaviour y = (3/2π)3/2 exp(−3n2/2). One sees that compared
to this reference the measured distributions are depleted for small n ≪ 1 (where the data does not
scale) and enhanced for n ≈ 1.
0.0 0.5 1.0 1.5 2.0 2.5
n=r/b
BFM s=31
BFM s=63
BFM s=127
BSM s=32
BSM s=64
BSM s=128
f(n) from Eq. (3)
Enhancement
epletion
FIG. 12: (Color online) Deviation δG(r, s) = G(r, s) − G0(r, s) of the measured segmental
size distribution from the Gaussian behavior G0(r, s) expected from Flory’s hypothesis for sev-
eral s and both models as indicated in the figure. As suggested by Eq. (3), we have plotted
y = (δG(r, s)/G0(r, s))/(ce/
s) as a function of n = r/be
s. The Gaussian reference distribution
has been computed according to Eq. (2) for the measured effective bond length be. A close to
perfect data collapse is found for both models. This shows that the deviation scales linearly with
u(s) ≈ ce/
s, as expected. The bold line indicates the universal function of f(n) predicted by
Eq. (3).
0.0 0.5 1.0 1.5 2.0 2.5
m=r/R
BFM s=31
BFM s=63
BFM s=127
BSM s=32
BSM s=64
BSM s=128
f(m) from Eq. (29)
Enhancement
letio
FIG. 13: (Color online) Replot of the relative deviation of the measured segment size distribution,
y = (δG(r, s)/G0(r, s)))/(ce/
s), as a function of m = r/Re(s). The figure highlights that the
measured segment size is the only length scale relevant for describing the deviation from Flory’s
hypothesis. The same data sets and symbols are used as in the previous Fig. 12.
v(k) v(k) v(k)~ ~ ~
q q q q q
q−k q−k
FIG. 14: Interaction diagrams used in reciprocal space for the calculation of δG(q, t) in the scale
free limit. There exist three nonzero contributions to first-order perturbation, the first involving
two points inside the segment (first two lines of Eq. (A4)), the second one point inside and one
outside the segment (third line of Eq. (A4)) and the third one point on either side of the segment
(last line of Eq. (A4)). Momentum q flows from one correlated point to the other. Integrals are
performed over the momentum k. Dotted lines denote the effective interactions ṽ(k) given by
Eq. (9), bold lines the propagators which carry each a momentum q or q − k as indicated.
[1] M. Rubinstein and R. Colby, Polymer Physics (Oxford University Press, 2003).
[2] P. J. Flory, J. Chem. Phys. 13, 453 (1945).
[3] P. J. Flory, J. Chem. Phys. 17, 303 (1949).
[4] P. J. Flory, Statistical Mechanics of Chain Molecules (Oxford University Press, New York,
1988).
[5] P. G. de Gennes, Scaling Concepts in Polymer Physics (Cornell University Press, Ithaca, New
York, 1979).
[6] M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Clarendon Press, 1986).
[7] S. Edwards, Proc. Phys. Soc. 85, 613 (1965).
[8] S. Edwards, Proc. Phys. Soc. 88, 265 (1966).
[9] S. Edwards, J.Phys.A: Math.Gen. 8, 1670 (1975).
[10] M. Muthukumar and S. Edwards, J. Chem. Phys. 76, 2720 (1982).
[11] L. Schäfer, Excluded Volume Effects in Polymer Solutions (Springer-Verlag, New York, 1999).
[12] A. N. Semenov and A. Johner, Eur. Phys. J. E 12, 469 (2003).
[13] A. N. Semenov and S. P. Obukhov, J. Phys.: Condens. Matter 17, 1747 (2005).
[14] P. Beckrich, A. Johner, A. N. Semenov, S. P. Obukhov, H. C. Benôıt, and J. P. Wittmer,
Macromolecules (2007), in press, cond-mat/0701261.
[15] P. Beckrich, Ph.D. thesis, Université Louis Pasteur, Strasbourg, France (2006).
[16] L. Schäfer, M. Müller, and K. Binder, Macromolecules 33, 4568 (2000).
[17] R. Auhl, R. Everaers, G. Grest, K. Kremer, and S. Plimpton, J. Chem. Phys. 119, 12718
(2003).
[18] J. P. Wittmer, H. Meyer, J. Baschnagel, A. Johner, S. P. Obukhov, L. Mattioni, M. Müller,
and A. N. Semenov, Phys. Rev. Lett. 93, 147801 (2004), cond-mat/0404457.
[19] J. P. Wittmer, P. Beckrich, A. Johner, A. N. Semenov, S. P. Obukhov, H. Meyer, and
J. Baschnagel, Europhys. Lett. 77, 56003 (2007), cond-mat/0611322.
[20] J. P. Wittmer, P. Beckrich, F. Crevel, C. C. Huang, A. Cavallo, T. Kreer, and H. Meyer,
Computer Physics Communications (2007), in press, cond-mat/0610359.
[21] A. Cavallo, M. Müller, J. P. Wittmer, and A. Johner, J. Phys.: Condens. Matter 17, 1697
(2005), cond-mat/0412373.
[22] H. Meyer, T. Kreer, A. Cavallo, J. P. Wittmer, and J. Baschnagel, Eur. Phys. J. Special Topics
141, 167 (2007), cond-mat/0609127.
[23] J. Baschnagel, J. P. Wittmer, and H. Meyer, in Computational Soft Matter: From Synthetic
Polymers to Proteins, edited by N. Attig (NIC Series, Jülich, 2004), vol. 23, pp. 83–140,
cond-mat/0407717.
[24] I. Carmesin and K. Kremer, Macromolecules 21, 2819 (1988).
[25] H. Deutsch and K. Binder, J. Chem. Phys 94, 2294 (1991).
[26] W. Paul, K. Binder, D. Heermann, and K. Kremer, J. Phys. II 1, 37 (1991).
[27] A. Kron, Polym. Sci. USSR 7, 1361 (1965).
[28] F. Wall and F. Mandel, J. Chem. Phys. 63, 4592 (1975).
[29] L. Mattioni, J. P. Wittmer, J. Baschnagel, J.-L. Barrat, and E. Luijten, Eur. Phys. J. E 10,
369 (2003), cond-mat/0212433.
[30] N. Karayiannis, A. Giannousaki, V. Mavrantzas, and D. Theodorou, J. Chem. Phys. 117,
5465 (2002).
[31] B. J. Banaszak and J. J. de Pablo, J. Chem. Phys. 119, 2456 (2003).
[32] M. Müller, J. P. Wittmer, and J.-L. Barrat, Europhys. Lett. 52, 406 (2000), cond-
mat/0006464.
[33] T. Kreer, J. Baschnagel, M. Müller, and K. Binder, Macromolecules 34, 1105 (2001).
[34] J. M. Deutsch, Phys. Rev. Lett. 54, 56 (1985).
[35] A. N. Semenov, in Theoretical Challenges in the Dynamics of Complex Fluids, edited by
T. McLeish (Kluwer, Dordrecht, 1997), pp. 63–86.
[36] H. Meyer and F. Müller-Plathe, J. Chem. Phys. 115, 7807 (2001).
[37] H. Meyer and F. Müller-Plathe, Macromolecules 35, 1241 (2002).
[38] K. Kremer and G. Grest, J. Chem. Phys. 92, 5057 (1990).
[39] M. Allen and D. Tildesley, Computer Simulation of Liquids (Oxford University Press, Oxford,
1994).
[40] O. Kratky and G. Porod, Rec. Trav. Chim. 68, 1106 (1949).
[41] J. Hansen and I. McDonald, Theory of simple liquids (Academic Press, New York, 1986).
[42] B. Alder and T. Wainwright, Phys. Rev. A 1, 18 (1970).
[43] P. Beckrich, J. P. Wittmer, H. Meyer, S. P. Obukhov, A. N. Semenov, and A. Johner (2007),
in preparation.
[44] A. N. Semenov, Journal de Physique II France 6, 1759 (1996).
[45] S. P. Obukhov and A. N. Semenov, Phys. Rev. Lett. 95, 038305 (2005).
[46] A. N. Semenov and M. Rubinstein, Eur. Phys. J. B 1, 87 (1998).
[47] J. Curro, K. S. Schweizer, G. S. Grest, and K. Kremer, J. Chem. Phys. 91, 1359 (1991).
[48] M. Fuchs, Z. Phys. B 5, 521 (1997).
[49] M. Schulz, H. Frisch, and P. Reineker, New Journal of Physics 6, 77/1 (2004).
[50] C. Likos, Physics Reports 348, 267 (2001).
[51] F. Eurich and P. Maass, J. Chem. Phys. 114, 7655 (2001).
[52] G. Yatsenko, E. J. Sambriski, M. A. Nemirovskaya, and M. Guenza, Phys. Rev. Lett. 93,
257803 (2004).
[53] J.-U. Sommer and K. Saalwächter, Eur. Phys. J. E 18, 167 (2005).
[54] C. Svaneborg, G. Grest, and R. Everaers, Europhys. Lett. 72, 760 (2005).
[55] R. Everaers, S. Sukumaran, G. Grest, C. Svaneborg, A. Sivasubramanian, and K. Kremer,
Science 303, 823 (2004).
[56] W. Paul and G. Smith, Rep. Progr. Phys. 67, 1117 (2004).
[57] P. Verdier, J. Chem. Phys. 45, 2118 (1966).
[58] R. Jones, S. Kumar, D. Ho, R. Briber, and T. Russell, Nature 400, 146 (1999).
[59] N. G. van Kampen, Stochastic processes in physics and chemistry (North-Holland, Amster-
dam, 1992).
[60] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions (Dover, New York,
1964).
[61] It is interesting to compare the numerical value I(∞) ≈ 1.59 obtained for the r.h.s of Eq. (16)
with the coefficients one would obtain by computing Eq. (15) either with the effective potential
ṽ(q) for infinite chains given by Eq. (9) or with the Padé approximation, Eq. (17). Within
these approximations of the full linear response formula, Eq. (14), the coefficients can be
obtained directly without numerical integration yielding overall similar values. In the first
case we obtain 15/8 ≈ 1.87 and in the second 11/8 ≈ 1.37. While the first value is clearly
not compatible with the measured end-to-end distances, the second yields a reasonable fit,
especially for small N < 1000, when the data is plotted as in Fig. 5. Ultimately, for very long
chains the correct coefficient should be 1.59 as is indicated by the dash-dotted line in Fig. 4.
[62] These topology non-conserving moves yield configurations which are not accessible with the
classical scheme with jumps in 6 directions only. Concerning the static properties we are
interested in this paper both system classes are practically equivalent. This has been confirmed
by comparing various static properties and by counting the number of monomers which become
“blocked” (in absolute space or with respect to an initial group of neighbor monomers) once
one returns to the original local scheme. Typically we find about 10 blocked monomers for a
system of 220 monomers. The relative difference of microstates is therefore tiny and irrelevant
for static properties. Care is needed, however, if the equilibrated configurations are used to
investigate the dynamics of the topology conserving BFM version. The same caveats arise for
the slithering snake and double-bridging moves.
[63] We have verified that the alternative definition 〈em=n+s · en〉 with en being the normalized
bond vector yields very similar results. This is due to the weak bond length fluctuations,
specifically at high densities, in both coarse-grained models under consideration. It is possi-
ble that other models show slightly different power law amplitudes cP depending on which
definition is taken.
	Flory's ideality hypothesis revisited
	Physical idea and sketch of the perturbation calculation
	Scaling arguments
	Perturbation calculation
	Computational models and technical details
	Bond fluctuation model
	Bead spring model
	Numerical results
	The swelling of chain segments
	Chain connectivity and recursion relation
	Intrachain bond-bond correlations 
	Higher moments and associated coefficients
	The segmental size distribution
	Conclusion 
	Acknowledgments
	Moments of the segmental size distribution and their generating function
	References
