Thermal Breakage and Self-Healing of a Polymer Chain under Tensile Stress

# Thermal Breakage and Self-Healing of a Polymer Chain under Tensile Stress

## Abstract

We consider the thermal breakage of a tethered polymer chain of discrete segments coupled by Morse potentials under constant tensile stress. The chain dynamics at the onset of fracture is studied analytically by Kramers-Langer multidimensional theory and by extensive Molecular Dynamics simulations in - and -space. Comparison with simulation data in one- and three dimensions demonstrates that the Kramers-Langer theory provides good qualitative description of the process of bond-scission as caused by a collective unstable mode. We derive distributions of the probability for scission over the successive bonds along the chain which reveal the influence of chain ends on rupture in good agreement with theory. The breakage time distribution of an individual bond is found to follow an exponential law as predicted by theory. Special attention is focused on the recombination (self-healing) of broken bonds. Theoretically derived expressions for the recombination time and distance distributions comply with MD observations and indicate that the energy barrier position crossing is not a good criterion for true rupture. It is shown that the fraction of self-healing bonds increases with rising temperature and friction.

###### pacs:
05.40.-a, 82.20.Uv, 82.20.Wt, 02.50.Ey,

## I Introduction

A great variety of problems both in material and basic science rely on the fundamental understanding of the intramolecular dynamics and kinetics of fragmentation (bond rupture) of linear macromolecular subject to a tensile force. Typical examples comprise material failure under stress Kausch (); Crist (), polymer rupture Garnier (); Gaub (); Saitta (); Maroja (); Rohrig (), adhesion Gersappe (), friction Klafter (), mechanochemistry Aktah (); Beyer (), and biological applications of dynamical force microscopy Harris (); Pereverzev (). In particular, the problem of polymer rupture as a kinetic process has a longstanding history and dates back to the publications of Bueche Bueche () and Zhurkov et al. Zhurkov (). In these papers the breaking of intramolecular bonds is treated as a thermally activated process and is described by Arrhenius’s formula for the rate of bond scission.

In recent years investigations have been complemented by ample application of computer experiments. Comprehensive Molecular Dynamic (MD) simulations of chain fragmentation (at constant strain) have been carried out, whereby harmonic Taylor (); Lee (), Morse Stember (); Bolton (); Sebastian (), or Lennard-Jones Oliveira_1 (); Oliveira_2 (); Oliveira_3 (); Oliveira_4 () chain models have been used. One of the principal questions to be answered is : “How long will it take for this system to break?”. A theoretical interpretation of MD-results Oliveira_3 (); Oliveira_4 (), based on an effectively one-particle model, has been suggested in terms of Kramers rate theory Kramers (). Thus, the activation energy of the one-particle model has been found to be close to the barrier height observed in the simulations, whereas the measured frequency of bond scission appears more than two orders of magnitudes smaller then the corresponding frequency, predicted by Kramers theory. The authors interpret this controversial result as a manifestation of the existence of collective modes which are missing in the one-particle Kramers theory.

It has also been noticed that a truly irreversible break may occur, if bonds are stretched to lengths, considerably larger than the one corresponding to the barrier position. This indicates the possibility for bond recombination whereby the chain integrity is restored with some finite probability.

The problem of polymer fragmentation has also been studied theoretically Sebastian () for the case of constant stress when a tethered chain of segments is subjected to a pulling force at the free chain end. The consideration has been based on a multidimensional version of the transition state theory (TST). Friction is then taken into account by coupling the polymer to a set of harmonic oscillators, simulating thus the presence of a thermostat. A comparison of the calculated breaking rate with the corresponding MD observation Oliveira_1 () shows again that the theoretical rate is about 250 times larger. Also the role of the bond healing process has been discussed which helps to improve to some extent the agreement between theory and simulation. Nevertheless, despite the multidimensional nature of TST, it does not take into account properly the collective unstable mode development, which leads, in our opinion, to the essential overestimation of the breaking rate.

The collectivity effect has been recently treated Sain () for constant strain and periodic boundary conditions (a ring polymer) on the basis of the multidimensional Kramers approach Langer (); Hanngi (). Within this approach the development of a collective unstable mode and the effect of dissipation can be described consistently. It has been shown that in this case the effective break frequency is of the same order of magnitude as the one observed in the simulation.

In the present work we develop this approach further for the case of a tethered Morse chain, consisting of segments and subjected to a constant tensile force applied at its free end. We derive analytic expressions for the scission rate of the bonds, its distribution along the polymer chain, and its variation with changing temperature and dissipation. For comparison with computer experiment, we also perform extensive MD simulations in both one- , and three dimensions, , and witness significant differences in the fragmentation behavior of the chain, despite the observed good agreement between theoretical predictions and simulation data. A major objective of the current study is the elucidation of the problem of bond recombination (self-healing) which has found little attention in literature so far. To this end we derive analytic expressions for the life times and extension distances of healing bonds and compare them to our MD results.

The paper is organized as follows. In Section II we give a sketch of the multidimensional Kramers-Langer escape theory Langer (); Hanngi (). We also outline the problem of multiple points of exit from potential well Matkowsky (); Gardiner () which is necessary for treating bond rupture with respect to the consecutive number of each bond along the chain. In Section III we present our model of one-dimensional Morse string of beads and consider the eigenvalue problem in the vicinity of the metastable minimum of the effective potential and at the barrier (saddle point) which is needed for the description of the unstable collective mode. Section IV gives briefly details of our MD simulation and presents the main numeric results as well as their interpretation in the light of our theoretical approach. In Section V the healing process is discussed in terms of distributions of healing times and bond extensions. We give also the theoretical interpretation of this process based on the solution of the Kramers equation Risken () using an inverted harmonic potential for representation of the barrier. We conclude in Section VI and outline some future developments.

## Ii Kramers-Langer multidimensional escape theory

### ii.1 Rate of escape

The calculation of the rate of escape is based on the -dimensional (in the total phase space , , Goldstein ()) Fokker-Planck equation Risken ()

 ∂P(x,t)∂t=2N∑i,j=1∂∂xiMij[∂H∂xjP(x,t)+T∂∂xjP(x,t)] (1)

for the probability distribution function . In eq. (1) the Hamiltonian has a general form . The -matrix , where is the matrix of Onsager coefficients and skew-symmetric matrix Goldstein ()

 →A=⎛⎝→0→1−→1→0⎞⎠ (2)

where is the unit matrix and is the zero matrix. The eq.(1) can be seen as a continuity equation, , where the probability current is given by .

It is assumed Langer () that there is a metastable state which is separated with a barrier from another stable state. The coordinates at the saddle point which separates these two states is denoted by . The escape from the metastable minima is a comparatively rare event, so that one can treat the process close to as a stationary one, i.e. . Thus one gets

 2N∑i,j=1∂∂xiMij[∑kESjk(xk−xSk)+T∂∂xj]P(x)=0 (3)

where also the harmonic approximation around has been used, i.e., the Hamiltonian reads

 H({xi})=ES+122N∑j,k=1ESjk(xj−xSj)(xk−xSk) (4)

where and the Hessian matrix at .

One should impose the following boundary conditions. Near the metastable state the distribution function is the equilibrium one, i.e.

 P(x)=Peq(x)=Z−1Aexp(−βH)atx≃{xAi} (5)

On the other hand, all states around the stable minimum (which is far beyond the top of the barrier!) are removed by a sink, i.e.,

 P(x)≃0at{xi}far beyond{xSi} (6)

With a transformation to new coordinates , which are principal-axis coordinates for the Hamiltonian given by eq.(4), one obtains

 H(ξ)=ES+122N∑n=1λnξ2n+… (7)

where the vector . The matrix is orthogonal one, i.e. and are the eigenvalues of the matrix . Since is a saddle point one of the ’s, say , is negative. The standard trick would be to look for the solution in the form , where is a new function. Thus in the -coordinate system the steady - state Fokker Planck equation for takes the form

 2N∑n,k=1(˜Γnk∂2W∂ξn∂ξk−βλnξn˜Mnk∂W∂ξk)=0 (8)

where .

In the same manner as for the one-dimensional Kramers problem Kramers () one could claim that the function depends only on a linear combination of all , i.e. , where . The prime in this expression indicates that we omit all ’s for which . The resulting equation reads: . As in the one-dimensional Kramers problem Kramers () the coefficient in front of is a linear function of , i.e. . As a result

 λn∑k˜MnkUk=κUn (9)

i.e., the coefficients are solutions of the eigenvalue problem eq. (9). Thus,

 Tλ+d2F(u)du2−udF(u)du=0 (10)

where

 λ+=1κ∑n,kUn˜ΓnkUk (11)

There is a simple physical interpretation of the eigenvalue problem given by eq. (9) Langer (). Namely, the equation of motion for the average value, can be written as

 ∂∂t⟨ξn(t)⟩=∫2N∏j=1dξj˜J({ξ},t)=−TZA∑k˜Mnk∫2N∏j=1dξj∂W∂ξke−βH (12)

With integration by parts and the harmonic approximation, eq. (7), one arrives at the following expression

 ∂∂t⟨ξn(t)⟩=−∑k˜Mnkλk⟨ξk(t)⟩ (13)

The unstable solution of this equation (which describes the decay of a metastable state) is given by a negative eigenvalue , namely

 ⟨ξn(t)⟩=Xne−κt (14)

Substitution of eq.(14) in eq. (13) leads to or

 λn∑k˜MnkλkXk=κλnXn (15)

This equation is identical to eigenvalue problem, eq. (9), provided that . There is a negative eigenvalue and as a result a negative (see eq.(11)) which corresponds to the unstable mode. This clear physical interpretation justifies the linear combination ansatz which has been used upon the derivation of eq. (10). We will show in Sec. IV devoted to MD-simulation results that the law given by eq. (14) actually holds for the breaking bonds.

With the negative eigenvalue, , the solution of eq.(10) takes the form

 F(u)=1√2π|λ+|T∫∞udzexp(−z22|λ+|T) (16)

In eq.(16) we take into account the boundary conditions, eqs.(5) and (6) which require (at the metastable well) and (around the stable minimum).

Consider the rate of the metastable state decay. The main quantity which should be used for this purpose is the the probability current

 ˜Jn=−TZA∑k˜MnkUkdF(u)due−βH (17)

To obtain the total probability flux over the barrier one should integrate the current, eq. (17), over the hypersurface containing the saddle point. The resulting flux reads

 J=∫2N∏i=1dξiδ(u)2N∑n=1Un˜Jn({ξi}) (18)

The calculation of this integral is given in Appendix of ref. Hanngi (). The calculation yields

 J=|κ|2π(2πT|λ1|)1/22N∏n=2(2πTλn)1/2e−βESZA (19)

In order to obtain the rate constant one must divide the flux over the popularion in the metastable well

 nA=∫Awell2N∏i=1dξiPeq({ξi})=e−βEAZA2N∏i=1(2πTλan)1/2 (20)

where we have used the harmonic approximation near the metastable state , i.e. . Combining eq.(19) and eq.(20), one finds for the rate constant, , the following result Langer ()

 k=|κ|2π[det(EA/2πT)|det(ES/2πT)|]1/2e−βEb (21)

where the activation barrier . We recall that are the eigenvalues (without zero-modes) of the Hessian matrix at the saddle point whereas are the the eigenvalues of the Hessian matrix at the metastable point .

### ii.2 The eigenvalue problem

Consider in more detail the eigenvalue problem given by eq. (15), i.e.. In the initial -space this equation reads

 ∑k,rMnkESkrXr=κXn (22)

where and the friction as well as the matrix are given as

 Γij=mγ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0⋯0⋮⋮0⋯00⋯0⋮⋱⋮0⋯00⋯0⋮⋱⋮0⋯01⋯0⋮⋱⋮0⋯1⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠Aij=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0⋯0⋮⋮0⋯01⋯0⋮⋱⋮0⋯1−1⋯0⋮⋱⋮0⋯−10⋯0⋮⋮0⋯0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ (23)

The Hessian matrix and the -dimensional column-vector are

 ESij=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝∂2V∂qi∂qj0⋯0⋮⋮0⋯00⋯0⋮⋮0⋯01/m⋯0⋮⋱⋮0⋯1/m⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠Xi=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝q1q2⋮qNp1p2⋮pN⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ (24)

After substitution of eqs. (23) and (24) into eq. (22), and exclusion of momentums , one arrives at the eigenvalue problem

 N∑j=1[VSij−mγκδij]N×N - matrixXj=−mκ2Xi (25)

where the notation is used. The corresponding characteristic equation Lancaster () reads

 det[mκ2δij−mγκδij+VSij]=N∏k=1[mκ2−mγκ+λk]=0 (26)

where is a set of eigenvalues of the Hessian . As mentioned before, only one eigenfunction, say , is negative. This eigenvalue determines then the equation for the transmission factor , i.e.

 κ2−γκ−|λ1|m=0 (27)

The negative solution reads

 κ=−√γ24+|λ1|m+γ4<0 (28)

Eq. (26) has been discussed first in ref. Weidenmuller ().

### ii.3 First-passage-time approach

The method of Section II.1 is a quite general approach to evaluate the rate of escape from the metastable state. The alternative to the flux-over-population method involves the concept of mean first-passage time. For an arbitrary stochastic multi-dimensional process the mean first-passage time (MFPT) is defined as the average time elapsed until the process starting out at point leaves a prescribed domain which includes the point of the reactant state. is sometimes called domain of attraction of the metastable state, and is shown in Fig. 1

For the FPT investigation it is convenient to use the Fokker-Planck equation (for the conditional probability of visiting the point at time , provided it starts from at ) which includes the so called adjoint operator and has the form Risken ()

 ∂∂tP(y,t;x,0)=L†(x)P(y,t;x,0) (29)

where the adjoint operator

 L†(x)=∑iKi(x)∂∂xi+ε∑i,jDij(x)∂2∂xi∂xj (30)

In the present case the drift and diffusion coefficients are respectively: and . In eq. (30) is responsible for the noise intensity (in our case plays the role of the temperature), so that at the noise is weak.

The -th order moments of the first-passage time, , may be iteratively expressed in terms of the adjoint operator, Risken ()

 L†(x)Tn(x) = −nTn−1(x),atx∈Ω Tn(x) = 0,atx∈∂Ω (31)

where the second equation simply means that the first passage time is zero, provided the trajectory starts at the separatrix . The hierarchy given by eq. (31) should be supplemented by the initial condition which is evident from the normalization condition for the FPT distribution. Therefore, the equation for the mean first passage time (MFPT), , reads

 L†(x)T1(x) = −1,atx∈Ω T1(x) = 0,atx∈∂Ω (32)

On the other hand, for small noise, , a trajectory starting within will typically first approach the attractor and stay within its neighborhood for a long time until an occasional fluctuation drives it to the separatrix . Hence, MFPT assumes the same large value everywhere in except for a thin layer along the boundary . Under these condition the solution of the hierarchical equation, eq.(31), can be obtained in the following form Talkner (): . In result the probability density of the first passage times reads Talkner ()

 p(t)=τ−1exp(−t/τ) (33)

(we recall that ). The validity of this PDF has been shown in our MD-simulation (see Sec. IV). Finally, it can be proven Talkner (); Hanngi () that the Kramers escape rate .

### ii.4 Distribution of exit points

Assume that there are a number of saddle points, where , (referred to as exit points)which lie on the separatrix. We study then the exit events (see Fig. 1 where only two exit points are shown). One may ask : “What is the probability to exit through a particular exit point irrespective of the time it takes?” This problem has been discussed first by Matkowsky and Schuss Matkowsky () and later by Gardiner Gardiner (). Here we give an explicit solution which may be compared to MD-simulations.

The probability for escape through an exit point if the random trajectory starts at is defined as

 π(bS;x)=∞∫0dt′∑iνi(bS)Ji(bS,t′|x,0) (34)

where is a component of the vector normal to the separatrix unit vector at pointing out of region . The flux counts only trajectories which start at in and approach the exit point at the boundary (separatrix) at time moment . The integral over time in eq. (34) implies that this probability is calculated irrespective of the time needed for escape.

The probability is governed by the backward stationary FPE, i.e. (see Sec. 5.4.2 in Gardiner ())

 L†(x)π(bS;x)=0 (35)

The boundary conditions are: and for any if , i.e.,

 π(bS;x)=δs(bS−x)forx∈∂Ω (36)

Using the same arguments as in Sec. II.3, one can show that in the small noise limit, , the probability remains the same everywhere inside apart from a thin layer along the boundary . Now we specify the general eq. (35) for the special case when the drift term has the form of potential (i.e., it is derivative of the Hamiltonian). Then eq. (35) reads

 ∑i,j[−Mij∂H∂xj∂∂xi+εMij∂2∂xi∂xj]π(bS;x)=0 (37)

Close to the saddle point the Hamiltonian can be treated in the harmonic approximation, eq. (4), and the drift velocity in eq. (37) becomes

 Ki=−∑jMij∂H∂xj=−∑j,kMijESjk(xk−bSk) (38)

Since the solution changes only within a thin layer along the boundary (in direction normal to the boundary layer!), one may introduce new local coordinates where measures the distance from and is the set of tangential variables measuring the orientation around . The coordinates and are chosen so that

 ∇z(u) = ν(u) ν(u)⋅∇yp(u) = 0 z(bS) = 0 (39)

where . The first equation in (39) means that the coordinate changes in direction of the -vector. The second equation implies that the coordinates are parallel to . In terms of new variables

 ∇iπ=νi∂π∂z+∑p∇iyp(x)∂π∂yp (40)

and

 ∇i∇jπ = νiνj∂2π∂z2+2∑pνi∇jyp(x)∂2π∂z∂yp+∑p,s∇iyp(x)∇jys(x)∂2π∂yp∂yr (41) + ∇i∇jz(x)∂π∂z+∑p∇i∇jyp(x)∂π∂yp

One should keep in mind that for the exit event at only one coordinate is relevant, namely, the one traversing the saddle point in direction of vector. This is our new -coordinate, which parameterizes the displacement from the saddle point as follows

 xj = bSj+zνj z = √ερ (42)

Substituting eqs.(40)- (42) in eqs. (37) and (38), and keeping only the lowest order in , leads to

 λ+∂2π∂ρ2−ρ∂π∂ρ=0 (43)

where

 λ+=1κ∑i,jνiΓijνj=mγκ. (44)

In (44)one has used and taken into account (Appendix A) as well as

 κ=∑i,j,kνiMijESjkνk (45)

The solution of eq.(43) has the form

 π(bS;u,ρ)=δ(u−bS)+[C∞−δ(u−bS)](2π|λ+|)1/2ρ∫0dzexp[−z22|λ+|] (46)

where we took into account the boundary condition (vector ) and at , i.e., well inside the region the solution is a constant as this should be for a weak noise. In order to fix the constant one may multiply eq. (37) by and after integrating over (using also the integration by parts and the Gauss theorem) derive the surface integral

 ∫∂ΩdS{−pst∑i,jνiMij∇jHπ+ε[pst∑i,jνiMij∇jπ−π∑i,jνiMij∇jpst]}=0 (47)

Given that , the 1-st and 3-rd terms in eq. (47) cancel each other and therefore

 ∫∂ΩdSpst∑i,jνiMij∇jπ=0 (48)

The gradient , calculated from eq.(46) at the boundary (i.e. at ), reads

 ∇jπ|ρ=0=νj[C∞−δ(u−bS)](2π|λ+|)1/2 (49)

Substitution of eq.(49) into eq. (48) leads finally to the result for the exit probability

 π(bS)=C∞=√κ(bS)e−H(bS)/ε∫∂Ω√κ(u)e−H(u)/ε=√κ(bS)e−H(bS)/εM∑r=1√κ(br)e−H(br)/ε (50)

In the last equality in eq. (50) the surface integral is replaced by a sum over all saddle (exit) points. This is possible because all eigenvalues of the Hessian are positive within the boundary hypersurface and hence has minimums at (where ). The exit probability given by eq. (50) will be compared in Section IV with our MD-simulation results.

Finally, we stress that the choice of the coordinate system given by eq. (39) so that only -direction is physically relevant is in complete agreement with the Kramers-Langer approach where the linear combination (see the paragraph after eq. (8)) plays the role of the relevant coordinate. Indeed, the first relationship from eq.(39) can be formally solved as

 z(x)=∑iνj(xj−bSj) (51)

which after transformation to the principal-axis coordinates, i.e., leads to (where ). Moreover, on the separatrix and one recovers the Kramers-Langer choice of a relevant coordinate (see the paragraph after eq. (8)) as linear combination of all .

## Iii Breakage of an one-dimensional string of beads

Here we describe chain breakage by means of the Kramers approach, sketched in Section II.

### iii.1 Model

We consider a tethered one-dimensional string of beads which experiences a tensile force at the free end as depicted in Fig.2. Successive beads are joined by bonds, governed by the Morse potential, , where and are parameters, measuring the bond strength and elasticity. The total potential energy is

 V({xi})=N∑i=1UM(xi−xi−1)−fxN (52)

where we set (see Fig. 2). Upon change of variables, , one gets

 V({xi})=N∑n=1[UM(yn)−fyn]=N∑n=1U(yn) (53)

so the combined one-bond potential then reads .

Direct analysis of the one-bond potential indicates that the positions of the (metastable) minimum and of the maximum are given by

 y−,+=1aln[21±√1−~f] (54)

where the dimensionless force . The activation energy (barrier height) is given by

 Eb=U(y+)−U(y−)=D⎧⎪ ⎪⎨⎪ ⎪⎩√1−~f+~f2ln[1−√1−~f1+√1−~f]⎫⎪ ⎪⎬⎪ ⎪⎭ (55)

One can easily verify that decreases with . Since the Kramers’ theory implies , the force should not be too large. The characteristic frequencies at the minimum and maximum of , and , are therefore

 Ω21,2=a2Dm[√1−~f±(1−~f)]. (56)

Fig. 3a illustrates the Morse potential as well as the modified one-bond Morse potential . In order to test the role of anharmonicity (recall that the Kramers-Langer theory uses only harmonic approximation) we have also tested in our MD-simulation the Double-Harmonic potential constructed piecewise as for and for and shown in Fig. 3b. The corresponding barrier heights dependences on the pulling force are shown in inserts.

### iii.2 Eigenvalues close to the metastable minimum

One can verify that the determinants of the Hessian matrix, entering eq.(21), may be calculated exactly for our one-dimensional model and even the total eigenvalue problem may be solved analytically.

The Hessian at the metastable minimum has a block-matrix form

 EA=∂2H∂xi∂xj∣∣∣A=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝VAij0⋯0⋮⋮0⋯00⋯0⋮⋮0⋯01/m⋯0⋮⋱⋮0⋯1/m⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ (57)

where is the -matrix of the potential energy second derivatives. It has the following tridiagonal structure

 VAij=mΩ21⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝2−10−12−10−12−1⋱⋱⋱⋱ 0−12−10−11⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ (58)

The correct diagonal and off-diagonal elements follow from the double derivation of . For example, the upper-left element whereas the bottom-right element reads . The only nonzero off-diagonal elements are .

The block-matrix structure makes it possible to calculate the determinant as Lancaster ()