{'fraction_non_alphanumeric': 0.03910519769326448, 'fraction_numerical': 0.01292827258678527, 'mean_word_length': 4.049291083832905, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics of stick-slip in peeling of an adhesive tape

We investigate the dynamics of peeling of an adhesive tape subjected to a constant pull speed. We derive the equations of motion for the angular speed of the roller tape, the peel angle and the pull force used in earlier investigations using a Lagrangian. Due to the constraint between the pull force, peel angle and the peel force, it falls into the category of differential-algebraic equations requiring an appropriate algorithm for its numerical solution. Using such a scheme, we show that stick-slip jumps emerge in a purely dynamical manner. Our detailed numerical study shows that these set of equations exhibit rich dynamics hitherto not reported. In particular, our analysis shows that inertia has considerable influence on the nature of the dynamics. Following studies in the Portevin-Le Chatelier effect, we suggest a phenomenological peel force function which includes the influence of the pull speed. This reproduces the decreasing nature of the rupture force with the pull speed observed in experiments. This rich dynamics is made transparent by using a set of approximations valid in different regimes of the parameter space. The approximate solutions capture major features of the exact numerical solutions and also produce reasonably accurate values for the various quantities of interest.

We investigate the dynamics of peeling of an adhesive tape subjected to a constant pull speed. We derive the equations of motion for the angular speed of the roller tape, the peel angle and the pull force used in earlier investigations using a Lagrangian. Due to the constraint between the pull force, peel angle and the peel force, it falls into the category of differential-algebraic equations requiring an appropriate algorithm for its numerical solution. Using such a scheme, we show that stick-slip jumps emerge in a purely dynamical manner. Our detailed numerical study shows that these set of equations exhibit rich dynamics hitherto not reported. In particular, our analysis shows that inertia has considerable influence on the nature of the dynamics. Following studies in the Portevin-Le Chatelier effect, we suggest a phenomenological peel force function which includes the influence of the pull speed. This reproduces the decreasing nature of the rupture force with the pull speed observed in experiments. This rich dynamics is made transparent by using a set of approximations valid in different regimes of the parameter space. The approximate solutions capture major features of the exact numerical solutions and also produce reasonably accurate values for the various quantities of interest.

I. INTRODUCTION
Peeling is a kind of fracture that has been studied experimentally in the context of adhesion and is a technologically important subject. Experimental studies on peeling of an adhesive tape mounted on a cylindrical roll are usually in constant pull speed condition [1,2,3,4,5,6]. More recently, constant load experiments have also been reported [3,7]. Early studies by Bikermann [5], Kaeble [6] have attempted to explain the results by considering the system as a fully elastic object. This is clearly inadequate as it ignores the viscoelastic nature of the glue at the contact surface and therefore cannot capture many important features of the dynamics. The first detailed experimental study of Maugis and Barquins [1] show stick-slip oscillations within a window of pull velocity with decreasing amplitude of the pull force as a function of the pull velocity. Further, these authors report that the pull force shows sinusoidal, sawtooth and highly irregular (chaotic as these authors refer to) wave patterns with increasing velocities. More recently, Gandur et al. have carried out a dynamical time series analysis of the force waveforms, as well as those of acoustic emission signals and report chaotic force waveforms at the upper end of the pull velocities [3]. One characteristic feature of the peeling process is that the experimental strain energy release rate shows two stable branches separated by an unstable branch. Stick-slip behavior is commonly observed in a number of systems such as jerky flow or the Portevin-Le Chatelier (PLC) effect [8], frictional sliding [9], and even earthquake dynamics is thought to result from stick-slip of tectonic plates [10]. Stick-slip is characterized by the system spending most part of the time in the stuck state and a short time in the slip state, and is usually seen in systems subjected to a constant response where the force developed in the system is measured by dynamically coupling the system to a measuring device. One common feature of such systems is that the force exhibits "negative flow rate characteristic" (NFRC). Models which attempt to explain the dynamics of such systems use the macroscopic phenomenological NFRC feature as an input, although the unstable region is not accessible. This is true for models dealing with the dynamics of the adhesive tape as well. To the best of our knowledge, there is no microscopic theory which predicts the origin of the NFRC macroscopic law except in the case of the PLC effect [11,12] (see below).
As there is a considerable similarity between the peeling of an adhesive tape and the PLC effect, it is useful to consider the similarities in some detail. The PLC effect refers to a type of plastic instability observed when samples of dilute alloys are deformed under constant cross head speeds [13]. The effect manifests itself in the form of a series of serrations in a range of applied strain rates and temperatures. This feature is much like the peeling of an adhesive tape. Other features common to these two situations are: abrupt onset of the large amplitude oscillations at low applied velocities with a gradually decreasing trend and NFRC, which in the PLC effect refers to the existence of negative strain rate sensitivity of the flow stress. In the case of the PLC effect, the physical origin of the negative strain rate sensitivity is attributed to the ageing of dislocations and their tearing away from the cloud of solute atoms. Recently, the origin of the negative SRS has been explicitly demonstrated as arising from competing time scales of pinning and unpinning in the Ananthakrishna's model [11,12]. In the case of adhesive tape, the origin of NFRC can be attributed to the viscoelastic behavior of the fluid. (Constant load and constant load rate experiments are possible in the PLC also.) While simple phenomenological models based on NFRC explain the generic features of the PLC effect [14], there appears to be some doubts if the equations of motion conventionally used in the present case of peeling are adequate to describe the velocity jumps [2,4]. Indeed, these equations of motion are singular and pose problems in the numerical solutions.
Apart from detailed experimental investigation of the peeling process, Maugis and Barquins [1], have also contributed substantially to the understanding of the dynamics of the peeling process. However, the first dynamical analysis is due to Hong and Yue [2] who use an "N" shaped function to mimic the dependence of the peel force on the rupture speed. They showed that the system of equations exhibits periodic and chaotic stick-slip oscillations. However, the jumps in the rupture speed are introduced externally once the rupture velocity exceeds the limit of stability [4,15]. Thus, the stick-slip oscillations are not obtained as a natural consequence of the equations of motion. Therefore, in our opinion the results presented in Ref. [2] are the artifacts of the numerical procedure followed. Ciccotti et al. [4] interpret the stick-slip jumps as catastrophes. Again, the belief that the jumps in the rupture velocity cannot be obtained from the equations of motion appears to be the motivation for introducing the action of discrete operators on the state of the system to interpret the stick-slip jumps [4], though they do not demonstrate the correctness of such a framework for the set of equations. Lastly, there are no reports that explain the decrease in the amplitude of the peel force with increasing pull speed as observed in experiments. As there is a general consensus that these equations of motion correctly describe the experimental system, a proper resolution of this question (on the absence of dynamical jumps in these equations) assumes importance.
The purpose of this paper is to show that the dynamics of stick-slip during peeling can be explained using a differential-algebraic scheme meant for such singular situations [16] and demonstrate the rich dynamics inherent to these equations. In what follows we first derive the equations of motion (used earlier [2]) by introducing an appropriate Lagrangian for the system. Then, we use an algorithm meant to solve differential-algebraic equations [16] and present the results of our simulations for various parameter values. One of our major findings is that inertia has a strong influence on the dynamics. In addition, following the dynamization scheme similar to the one used in the context of the PLC effect [14], we suggest that the peel force depends on the applied velocity. Using this form of peel force leads to the decreasing nature of the magnitude of the pull force as a function of applied velocity. For certain values of the inertia, we find canard type solutions. These numerical results are captured to a reasonable accuracy using a set of approximations valid in different regimes of the parameter space. Even though, our emphasis is on demonstrating the correctness of these equations of motion and richness of the inherent dynamics that capture the qualitative features of the peeling process, we also attempt to make a comparison of the experimental results mentioned above to the extent possible.

II. EQUATIONS OF MOTION
For the sake of completeness, we start by considering the geometry of the experimental setup shown schematically in Fig. 1. An adhesive roll of radius R is mounted on an axis passing through O normal to the paper and is pulled at a constant velocity V by a motor positioned at O ′ with a force F acting along P O ′ . Let the distance between O and O ′ be l, and that between the contact point P to O ′ be L. The point P moves with a local velocity v which can undergo rapid bursts in the velocity during rupture. The force required to peel the tape is usually called the force of adhesion denoted by f . The two measured branches referred to earlier, are those of the function f in a steady state situation of constant pulling velocity (i.e., there are no accelerations). The line L makes an angle θ with the tangent at the contact point P . The point P subtends an angle α at O, with the horizontal line OO ′ . We denote the elastic constant of the adhesive tape by k, the elastic displacement of the tape by u, the angular velocity by ω and the moment of inertia of the roll by I. The angular velocity itself is identified by ω =α + v/R. The geometry of the setup gives L cos θ = −l sin α and L sin θ = l cos α − R which further gives, L 2 = l 2 +R 2 −2lR cos α. The total velocity V at O ′ is then made up of three contributions [1], given Following standard methods in mechanics, it is straightforward to derive the equations of motion for α and ω by considering (α,α, u,u) as the generalized co-ordinates. The corresponding Lagrangian of the system can be written as We write the dissipation function as where f (v, V ) physically represents the peel force which we assume is dependent on rupture speed as well as the pull speed assumed to be derivable from a potential function Φ(v, V ). The physical origin of this is due to the competition between the internal relaxation time scale of the viscoelastic fluid and the time scale determined by the applied velocity [17]. When the applied velocity is low, there is sufficient time for the viscoelastic fluid to relax. As we increase the applied velocity, the relaxation of the fluid gets increasingly difficult and thus behaves much like an elastic substance. The effect of competing time scales is well represented by Deborah number [18] which is the ratio of time scale for structural relaxation to the characteristic time scale for deformation. Indeed, in the studies on Hele-Shaw cell with mud as the viscous fluid, one observes a transition from viscous fingering to viscoelastic fracturing [19] with increasing rate of invasion of the displacing fluid.
As stated in the Introduction, the existing models do not explain the decreasing amplitude of pull force. Similar feature observed in the PLC serrations has been modeled using a scheme referred to as dynamization of the negative strain rate sensitivity (SRS) of the flow stress f (ǫ P ) [14,20], whereǫ p is the plastic strain rate. Based on arguments similar to the preceding paragraph, they modify this function to depend on the applied strain rate, ǫ a , i.e., the negative SRS of the flow stress is taken to be f (ǫ P ,ǫ a ) such that the gap between the maximum and the minimum of the function f (ǫ p ,ǫ a ) decreases with increasingǫ a . Following this, we consider f to depend on V also, in a way that the gap in f decreases as a function of the pull speed V (Fig. 2).
Using the Lagrange equations of motion, d dt we obtain the same set of ordinary differential equations as in Ref. [2] given bẏ with an algebraic constraint (The last equation results from the elimination of two second order equations for α.) In Eqs. (7), (9), and (10) we have used cosθ ≃ −sinα ∼ −α. While Eqs. (6)-(9) are differential equations, Eq. (10) is an algebraic constraint necessitating the use of differential-algebraic scheme to obtain the numerical solution [16]. The fixed point of Eqs. (6), (7), (9), and (10) is given by α = 0, ω = V /R, v = V, F = f (V, V ). (For numerical solution, in the above equations we have actually used sinα in place of α.) This point is stable for f ′ (V, V ) > 0 and unstable for f ′ (V, V ) < 0. As V is varied such that the sign of f ′ (V, V ) changes from negative to positive value, the system undergoes a Hopf bifurcation and a limit cycle appears. The limit cycles reflect the abrupt jumps between the two positive slope branches of the function f (v, V ).

III. ALGORITHM
The singular nature of these equations becomes clear if one were to consider the differential form of Eq. (10) where the prime denotes the derivative with respect to v. Equations (11) with Eq. (6), (7), and (8) [or (9)] constitute the full set of evolution equations for the vector (α, ω, F, v). However, it is clearly singular at points of extremum of f (v, V ), requiring an appropriate numerical algorithm. We note that Eqs.(6), (7), (8), and (10) can be written as where X = (α, ω, F, v), φ is a vector function that governs the evolution of X and M is a singular "mass matrix" [16] given by, Equation (13) is a differential-algebraic equation (DAE) and can be solved using the so called singular perturbation technique [16] in which the singular matrix M is perturbed by adding a small constant ǫ such that the singularity is removed. The resulting equations can then be solved numerically and the limit solution obtained as ǫ → 0. We have checked the numerical solutions for ǫ values ranging from 10 −7 to 10 −15 in some cases and the results do not depend on the value of ǫ used as long as it is small. The results presented below, however, are for ǫ = 10 −7 . We have solved Eq. (13) using a standard variable-order solver, MATLAB ODE15S program.
to give values of the extremum of the peel velocity that mimic the general form of the experimental curves [1]. The measured strain energy release rate G(V ) from stationary state measurements is shown in Fig. 2(b). The decreasing nature of the gap between the maximum and minimum of f (v, V ) for increasing V is clear from Fig.  2(a).
[The values of f (v, V ) could not be correctly determined as G(V ) is in J/m 2 requiring more details. However, the value of F max is closer to Ref. [2] and the jumps in v are similar to those in experiments.] The reason for using the form given by Eq. (14) is that the effects of dynamization are easily included through its dependence on the pulling velocity while more complicated terms are required to mimic completely the experimental curve (particularly the flat portion). However, we stress that the trend of the results remains unaffected when the actual experimental curve is used except for the magnitude of velocity jumps and the force values.

IV. RESULTS
We have studied the dynamics of the system of equations for a wide range of values of the parameters. We have found that transients for some regions of parameters space take considerable time to die out. The results reported here are obtained after these long transients are omitted. These equations exhibit rich dynamics, some even unanticipated. Here we report typical results for two important parameters, namely the pull velocity V (m/s) and the inertia I (kg m 2 ), keeping the elastic constant of the tape k = 1000 N/m, R = 0.1 m and l = 1 m [2]. The influence of k will also be mentioned briefly. ( Henceforth, we drop the units for the sake of brevity.) We find that the observed jumps of the orbit in the v-F plane occur in a fully dynamical way. More importantly, we find all the three possibilities namely, the orbit can jump when it approaches the limit of stability, before or beyond that permitted by f (v, V ). The dynamics can be broadly classified into low, intermediate and high regimes of inertia.
(i) Low inertia. Here also, there are three regimes: low, intermediate, and high pull velocity.
(a) Consider keeping inertia I at a low value (say I = 10 −5 ) and V also at a low value (say, near the top, say V = 0.4 ). Here we observe regular saw tooth form for the pull force F . The phase plot in the F -v plane is as shown in Fig. 3(a). The corresponding function f (v, V ) is also shown by the continuous curve. We see that the trajectory jumps almost instantaneously from B to C on reaching the maximum of f (v, V ) (or from D to A when it reaches the minimum). The system spends considerably more time on AB compared to that on CD. However, this feature of jumping of the trajectory at the limit of stability is only true for small values of I and when V is near the limit of stability. At slightly higher pull velocity, say V = 1, even for small I, say I = 10 −5 , the jumps occur even before reaching the top or bottom ( the points B and D) as can be seen from Fig. 3(b) for V = 1. The small amplitude high frequency oscillations seen in the phase plots [ Fig. 3(a), and 3(b)] on the branch AB are due to the inertial effect, i.e., finite value of I. These oscillations are better seen on the α(t) plot shown in Fig. 3(c). For these values of parameters, the system is aperiodic. b) As we increase V , even as the saw tooth form of F is retained, various types of periodic orbits [period 4 shown in Fig. 3(d) for V = 2] as well as irregular orbits are seen. In both cases (periodic as well as chaotic) the trajectory jumps from high velocity branch (CD) to the low velocity branch before traversing the entire branch or sometimes going beyond the values permitted by f . The value of F at which the orbit jumps is different for different cycles. For I = 10 −5 , at high velocity, say V = 4, the phase plot is periodic.
(ii) Intermediate and high inertia.  than f max (f min ).] When the trajectory jumps from AB to CD at the highest value of F for the cycle, the trajectory stays on CD for a significantly shorter time compared to the small inertia case (I = 10 −5 ) and jumps back to AB well before F has reached the minimum of f (v, V ), i.e., ∆F is much smaller than f max − f min . The pull force F cascades down through a series of back and forth jumps between the two branches till the lowest value of F for the cycle is reached. Note that F at the point n is less than f min . For the sake of clarity, two different portions of the trajectory are marked abcdef g and ijklmna corresponding to the top and bottom regions of the plot. The corresponding points are also identified on the F (t) plot. After reaching n, the orbit jumps to a on AB, the trajectory decides to move up all the way till F reaches a maximum value (larger than f max , the point b) without jumping to the CD branch. This part of F as a function of time, which is nearly linear on AB, (i.e., the segment ab) displays a noticeable sinusoidal modulation. The sinusoidal form is better seen in α [ Fig. 4(b)]. Note that the successive drops in F are of increasing magnitude. The jumps between the two branches in the v − F plane are seen as bursts of v [ Fig. 4(d)]. For these values of parameters, the system is periodic.
(b) As we increase V , the sinusoidal nature of F and α becomes more clear with its range becoming larger reaching a nearly sinusoidal at V = 4 for large I. [The range ab in Fig. 4(c) expands. Compare Fig. 5(a).] The magnitude of ∆F on the CD branch for small V and moderately or large I, gradually decreases with increasing V . The magnitude of ∆F itself decreases as I is increased. In the limit of large V and I, the drops in F and α become quite small which are now located near the maxima and minima of these curves. This is shown in Figs. 5(a) and 5(b). The sinusoidal nature is now obvious even in F (t) unlike for smaller V and I where it is clear only in α(t) for the low v branch. Note that for V = 4, the nature of f (v, 4) is nearly flat. This induces certain changes in the v − F phase plot that are not apparent in F and α. The jumps between the two branches are now concentrated in a dense band at low and high values of F . In this case, the maximum (minimum) value of F is significantly larger (less) than f max (f min ). These rapid jumps between the branches manifest as jitter at the top and bottom of F and α.
Unlike for small V  df /dv = 0, even when the maximum (minimum) of F are higher (lower) than that allowed by the stationary curve f (v, V ). The variables are aperiodic for the set of parameters. The phase plots appear to be generated by an effective f (v, V ) that is being cycled. [This visual feeling is mainly due to the fact that jumps between the branches still occur at the maximum and minimum of The influence of k is generally to increase the range of the pull force F as can be easily anticipated and to decrease the associated time scale.
It may be desirable to comment on the similarity of the nature of the force waveforms displayed by the model equations with those seen in experiments. As mentioned in the introduction, apart from qualitative statements on the waveforms in Ref. [1] (such as periodic, sawtooth etc., which are seen in the model as well), it should be stressed that there is a paucity of quantitative characterization of the waveforms. In this respect, the study by Gandur et al. [3] fills the gap to some extent. These authors have carried out a dynamical analysis of the time series for various values of the pull velocities (for a fixed value of the inertia corresponding to their experimental roller tape geometry). In order to compare this result, we have calculated the largest Lyapunov exponent for a range of values of I and V . The region of chaos is in the domain of small pull velocities V when I is small. The maximum Lyapunov exponent turns out to be rather high, typically around 7.5 bits/s in contrast to the small values reported in Ref. [2]. The large magnitude of the positive exponent in our case can be traced to the large changes in the Jacobian, as df (v, V )/dv varies over several order of magnitude(∼ 10 6 ) as a function of the peeling velocity and hence as a function of time. In contrast, Hong et al. use an N shaped curve where df (v, V )/dv is constant (and small) on both low and high V branches. However, these large values of Lyapunov exponents are consistent with rather high values obtained by Gandur et al. [3] from time series analysis of the pull force. We also find chaos for intermediate and high inertia in the region of high velocities where the value of the Lyapunov exponent is small, typically 0.5. The small value here again can be traced to the small changes in df (v, V )/dv at high velocities.
It must be mentioned that comparison with experiments is further complicated due to the presence of a two parameter family of solutions strongly dependent on both I and V . Thus, the phase diagram is complicated, i.e., the sequence of solutions encountered in the I-V plane as we change V or I or both does not in general display any specific ordering of periodic and chaotic trajectories (see Fig. 1 of Ref. [21]) usually found in the well known routes to chaos. (For instance 2 n periods should be observed before the odd periods [22].) Indeed, in our model, we find the odd periods 3,5,7 etc, on increasing V (or I), without seeing all the 2 n periods. (These odd periods also imply chaos at parameter values prior to that corresponding to these periods.) In view of this, a correct One quantitative result that can be compared with experiment is the decreasing trend of the force drop magnitude. We have calculated the magnitude of the force drops during stick-slip phase as a function the pull velocity V for both low ( I = 10 −5 ) and high (I = 10 −2 ) inertia cases. Figure 6 shows the monotonically decreasing trend of average ∆F (t) as V is increased, for both small and large I, a feature observed in experiments [1]. These two distinct behaviors are a result of the dynamization of f (v, V ) as in Eq. (14).
Finally, as another illustration of the richness of the dynamics seen in our numerical simulations, we show in Fig. 7, a plot of an orbit that sticks to unstable part of the manifold before jumping back to the AB branch. Such solutions are known as canards [23]. Though canard type of solutions are rare, we have observed them for high values of I and low values of V . In our case, such solutions are due to the competition of time scale due to inertia and that due to v. This again illustrates the influence of inertia of the roll on the dynamics of peeling.
It is clear that these equations exhibit rich and complex dynamics. A few of these features are easily understandable, but others are not. For instance, the sawtooth form of F for low inertia and low pull velocity can be explained as resulting from the trajectory sticking to stable part of f (v, V ) and jumping only when it reaches the limit of stability. For these parameter values, as the time spent by the system is negligible during the jumps between the branches AB and CD (and vice versa), the system spends most of the time on the branch AB and much less on CD due to its steep nature. Then, from Eq. (9), it is clear that we should find a sawtooth form whenever the peel velocity v jumps across the branch to a value of v larger than the pull velocity V .
However, several features exhibited by these system of equations are much too complicated to understand. We first list the issues that need to be explained.
(I) Small I. (a) We find high frequency tiny oscillations superposed on the linearly increasing F [on the AB branch or better seen in the α plot Fig. 3(c)]. This needs to be understood.
(b) The numerical solutions show that the influence of inertia can be important even for small I and small V . For instance, the jumps between AB and CD branches occur even before F reaches the extremum values of f .
(II) For intermediate and high values of inertia, for low V case. (a) We observe several relatively small amplitude saw tooth form of F on the descending part of the pull force F . These appear as a sequence of jumps between the two branches in the v − F plane which we shall refer to as the "jumping mode". A proper estimate of the magnitude of ∆F is desirable.
(b) In addition, there appears to be a critical value of F for a given cycle below which the return jumps from AB to CD stop and one observes a monotonically increasing trend in F [ab in Fig. 4(c)].
(III) High I and high V . (a) The jumps between the branches occur at a very high frequency [ Fig. 5(c)] and now are located near the extremum values of F and α. But these regions are separated by a stretch where the orbit monotonically increases on the AB branch and monotonically decreases on the CD branch. We need to elucidate the underlying causes leading to the switching between the jumping mode and monotonically increasing or decreasing mode.
(b) For large V , say V = 4 and large I (Fig. 5), the extent of values of F (t) range between 185 and 450 much beyond f (v, 4) whose range is around 300. This feature is less dominant for small I and small V case.

V. APPROXIMATE ANALYSIS OF THE DYNAMICS
As the dynamics is described by a coupled set of differential equations with an algebraic constraint, the results are not transparent. We first attempt to get insight into the complex dynamics through some simple approximations valid in each of the regimes of the parameters. Solution of these approximate equations will require appropriate initial values for the relevant variables which will be provided from the exact numerical solutions. Due to the nature of approximations, the results are expected to capture only the trend and order of magnitudes of the effects that are being calculated. But as we will show, even the numbers obtained match quite closely with the exact numerical results.
Our idea is to capture the dynamics through a single equation ( as far as possible or at most two as in the high I and V case) by including all the relevant time scales and solve the relevant equation on each branch. For this we note that the equation for α and v play a crucial role as the inertial contribution appears only through Eqs. (6) and (7) and the time spent by the system is controlled by the equation for v, Eq. (12). Using Eqs. (6) and (7), we getα The general equation for α can be written down by using Eq. (12), in Eq. (15), we geẗ In obtaining Eq. (17), we have used 1 + α ≃ 1 which is valid except for high I and high V . Further, in most cases, we can drop Rαα as the magnitude of this term is small and useḞ ≃ k(V − v). To be consistent we use We note however that even for high I and high V where α is not small, dropping 1+α and Rαα causes only 10% error.
Case I, small I On the low velocity branch AB, as v/R is small in Eq. (6), we can dropv term in Eq. 15. Thus, Note that for the low inertia case, sinα ≈ α approximation is clearly justified [see Eq. (7)]. Using this equation, we first get an idea of the relevant time scales as I is increased.
Case a Consider the low velocity branch AB where the small amplitude high frequency oscillations are seen on the nearly linearly increasing part of F [given by F (t) = F min + k(V − v)t, see for instance Fig. 3(b)]. A rough estimate of this time spent on this branch is obtained by (f max − f min )/kV ∼ t. Using f max ∼ 284 and f min ∼ 200, [from Fig. 3(b)], we get t = 0.084 (compared to the correct value of 0.063 which we shall obtain soon) which is much larger than the period of the high frequency oscillation. Thus, we could take the local value F for the purpose of calculating the period of the high frequency oscillation. Consider the orbit at the lowest value of F for which we can use F min ∼ f min (v, 1) ∼ 200. Then using Eq. (18), the frequency ν = (F R/I)/2π = 225 for I = 10 −5 which gives the period of oscillation T = 4.44 × 10 −3 . This agrees very well with the exact numerical value T = 4.1 × 10 −3 . This frequency decreases when the force reaches the maximum value F max ∼ f max (v, 1) ∼ 284 to ν = 261 giving T = 3.69 × 10 −3 which is again surprisingly close to the numerical value 3.72 × 10 −3 . In the numerical solutions, we find that the period gradually decreases [see Fig. 3(c)]. This feature is also easily recovered by using F = F min + k(V − v)t. This leads to an additional term in the equation of motion for α in Eq. (18), (19) where t is the time required for F to reach F max starting from F min . Here again the v term can be dropped. If F min was absent, the equation has the Airy's form. (Note that for this case also we could assume F min ∼ f min and F max ∼ f max .) Though this equation does not have an exact solution, we note that we could take α to have a sinusoidal form with 2πν = (F R/I) where F is treated as a slowly increasing parameter. (This assumption works quite well.) The above equation captures the essential features of the numerical solution. The numerical solution of Eq. (19) ( as also this representation ) gives the decreasing trend of the small amplitude high frequency oscillations. (Note that the Airy equation itself gives a decreasing amplitude [24]. ) We note that Eq. (18) is valid on the AB branch where v is small even for high inertia and small V case. Thus, we may be able to recover the gross time scales using this equation. Our numerical results show that as we increase the inertia, α exhibits a sinusoidal form on the AB branch [see Fig. 4(b)], although one full cycle is not seen. We note that though the value of α is much larger than that for small I, we can still use the above equations [Eq. (18) and (19)]. On this branch F increases from a value F min ∼ f min to a maximum F max ∼ f max . For large I = 10 −2 ( and V = 1), we get a rough estimate of the period by using the mean value of F ∼ 240 in Eq. 18. This gives a period T = 0.128 which already agrees satisfactorily with the numerically exact value T =0.11 considering the approximation used (i.e., using the mean F ). A better estimate can be obtained by using Eq. (19).
For the high I and V case, Fig. 5 for V = 4 shows that the wave forms are nearly sinusoidal except for a jitter at the top and bottom. For this case, f (v, 4) is nearly flat over the entire range of values of v, with a value ∼ 300. Here, even on the AB branch, we can not ignore thev term in Eq. (15). However, one sees that as v min = 0.335 and v max = 1.25 which suggest that to the leading order, we could ignore thev term. This gives the period T = 0.115. From Fig. 5(b), considering only the monotonically decreasing part (ab), the value of T /2 = 0.53 read off from the figure compares reasonably well with this value.
For the CD branch, as v is not small, thev/R term appears to be important in Eq. (15). Some idea of when this term is important can be had by looking at the time scales arising from inertia, namely, F R/I and the coefficient of the damping term, F/Rf ′ in Eq. (17). Consider V = 1 for I = 10 −5 and 10 −2 . The period obtained by assuming the mean value of F = 240 in F R/I gives 4 × 10 −3 for I = 10 −5 compared to 0.128 for I = 10 −2 . These numbers can be compared with the time scale Rf ′ /F which is 0.01 (where we have used f ′ ∼ 25 from numerical simulations for V = 1). This shows that for high inertia the damping coefficient F/Rf ′ in Eq. (17) is important. We will discuss this issue in more detail later.

Case b
Now we focus on the origin of jumps between the branches. We note that the jumps from CD to AB (or vice versa) occur only when the peel velocity v reaches a value where f ′ (v, V ) = 0. This also means that the time scale on each branch, whether it spends only a short time or not, is controlled by the equation for v. However, clearly the influence of inertia needs to be included. Here we present an approximate equation for v which is valid in the various limits of the parameters: where the time t is time spent on the branch considered ( low or high V ). In Eq. (21), we have again useḋ t with the same approximation used in Eq. (17). We now attempt to obtain correct estimates of the time spent by the orbit on each branch starting with the least complicated situation of the low inertia and small V . For this case, on the low velocity branch, one can use the sinusoidal solution for α, namely α = α in sin(2πνt + φ), where φ is a phase factor which also includes the contribution arising from the jump as well and 2πν = (F R/I) with F ≃ F in + kV t . Both α in and φ needs to be supplied. Alternately, one can use Eq. (18) with Eq. (21) for which we provide α in andα in at the point from the exact numerical solutions. We stress that this procedure is not equivalent to solving all the equations, as the only equation we use is Eq. (21) with the form of α already determined from the equation for α. [We note here that though we have used the sinusoidal form of α along with the initial conditions on α in , φ, it is simpler to supply the initial conditions α in ,α in and use Eq. (18). ] We note here that f ′ is a crucial factor that determines the time at which the orbit jumps from one branch to the other. Equation (21) needs to be integrated from v in to v f that are determined by the pulling velocity V , i.e., the form of f (v, V ).
For the low v branch f ′ term makes a significant contribution for the time spent by the trajectory on AB. Indeed, one can obtain the order of magnitude of the time spent by the orbit on AB by using a crude approximation for f ′ (v, 1) = −230(0.5 − v). This can be easily  Fig. 8. Using ∆t in F = F in + k(V − v)∆t gives ∆F = 63 and F = 274.5 which is in good agreement with the exact numerical value of F max = 275. It is interesting to note that this value is much less than f max = 283 [see Fig. 3(b)] or equivalently ∆F is less than f max − f min , what is also observed in our exact numerical simulation. The underlying mechanism of jumping of the orbit before F reaches f max also becomes clear from the analysis (Fig. 8). We note that the magnitude of the oscillatory component in v grows till it reaches v max permitted by f (v, 1). Then, the orbit has to jump to CD. Thus, the approximate solution gives an insight into the cause of the orbit jumping even before F reaches f max (for small I).
For the CD branch also, the dominant term is f ′ . Indeed, any reasonable function which has the same geometrical form of f shown in Fig. 2 will give good results for ∆t. Using the correct form of f ′ , we get ∆t = 0.005 which is close to the exact result. This again gives correct magnitude of ∆F = 72.5. In addition the nature of the v(t) obtained by this approximation is close to the exact numerical solution shown in the inset of Fig. 8.
Case II, intermediate and high I and low V The most difficult feature of our numerical solutions to understand is the dynamical mechanism leading to a series of drops in the pull force seen on the descending branch of F (t) for intermediate and high values of inertia and for a range of V values. Consider the high inertia and low V case (say I = 10 −2 and V = 1) shown in Fig. 4. As stated earlier, there are two different issues that need to be understood here. First, the series of small force drops ∆F and second the monotonic increasing nature of F on the AB branch.
In this case, as already discussed, the coefficient ofα, namely, F/Rf ′ term in Eq. (17) determines the time scale on CD, while on AB, the term F Rα/I dominates. Thus, the general equation valid for this case is where we use [Note that we have dropped k(V − v) term from Eq. (17) as this term does not have any dependence on α orα.] We start with the cascading effect. Consider the orbit when it is at the highest value of F in = 295.6 on the CD branch for which we can drop F Rα/I term. As f ′ is a function of v, and F also depends on time, it appears that we need to use coupled equationsα = −F α/Rf ′ with Eq. (21). However, the numerical solution of these equations show that one can make further approximation by taking f ′ to be constant taken at v = 15.54 and F = F in , as the time spent on this branch is very small. The error in using this approximation is within 10%. Indeed, using α in = −0.0304,α in = −160 and numerically integrating Eq. (21), along with Eq. (22) from v in = 15.54 to v f = 8.7 gives ∆t = 1.78 × 10 −3 . This compares reasonably with the numerical value of 1.89 × 10 −3 . Using this we get ∆F = 19.4 which compares well with the numerically exact value 19.8. At this point the orbit jumps to the low velocity AB branch (to the point e). Thus, as ∆t is small, for all practical purposes, we can ignore the dependence of f on v and F on t and use α to be an exponentially decreasing function for analytical estimates. These analytical estimates already give reasonably accurate numbers.
On the AB branch, the dominant time scale is determined by F Rα/I, and we can use the approximate sinusoidal form in Eq. (21), or Eq. (18) along with Eq. (21) for the time evolution from the point e. Integrating from v in = 0.0188 to v f = 0.4 with the appropriate initial values α in = 0.239,α in = 11.9 (or α, φ) and F in = 276.53, gives ∆t = 0.016 which again compares very well with exact numerical value ∆t = 0.0164. This gives ∆F = 11.61. The procedure for calculating the time spent by the orbit on CD and AB is the same and we find that successive values of ∆F increases which is again consistent with what is seen in Figs. 4(a) and 4(c).
Continuing this procedure, we find that a minimum value of F = 186.95 for the cycle is reached. Now consider the time evolution of F on AB that should lead to a monotonically increasing nature as seen in the numerically exact solution. As this point corresponds to the point at which the dynamics switches from the jumping mode to the monotonically increasing nature of F (i.e., the stretch ab), we discuss this in some detail. For the point a, we have used the initial condition α in = 0.0599,α in = 9.7 and integrating Eq. (21) and Eq. (18) (or the sinusoidal form of α) from v = v min = 0.0188 to v = v max = 0.4 gives ∆t = 0.117. This is nearly the value 0.114 obtained from the exact numerical integration. This gives ∆F = 117 and F max = 303.95 which compares very well with the exact numerical value. In addition, the growth form of v obtained from this approximation (continuous line) agrees very well with that of the exact numerical solution (dotted line) as shown in Fig. 9. The discrepancy seen in the figure can be reduced for instance if we include the terms neglected in Eq. (16) such asḞ and using 1 + α in Eq. (20). Now we come to the crucial question. How does the system know that it has to go from a to b, while just during the previous visit to the point k on AB branch lead only to a small increase in ∆F [ Fig. 4(a) and 4(c)] before jumping to CD? To understand this, we recall that on AB, a sinusoidal solution is allowed. First, one can notice a few differences in the initial conditions between the point a and k. For the point k, the initial conditions taken from the exact numerical solution are α in = 0.298 andα in = 18.3. ( F in = 193.37), while for the point a, α in = 0.0589,α in = 9.7. However, for α to begin a sinusoidal form, the initial value ofα = 18.3 is much higher than the natural slope. The local slope for any sinusoidal form is maximum when the variable is close to zero. In Fig. 4(b), the sinusoidal form starts when α is close to zero (∼ 0.0589 at a) where the local slope should be close to the maximum value. Near α ∼ 0, the local slope is the product of the maximum amplitude of α, say, α 0 (in the sinusoidal stretch ab) and 2πν. (We have assumed α = α 0 sin2πνt by dropping the phase factor.) Thus, one should haveα ≃ 2πνα 0 when α ∼ 0. Using the value α 0 = 0.23 from exact numerical solution and ν ∼ 6.88 at F in = 186.95, we find thaṫ α ∼ 10 near α ∼ 0. Indeed, this is satisfied only at a whereα in = 9.7. [Note that α is not symmetric around zero due to the presence of v in Eq. (6) which has been ignored for the purpose of present discussion.] However, α in = 18.3 at k is significantly higher than the slope permitted for α to start a sinusoidal sojourn. This forces the orbit to make one more small loop (AB to CD and back) so that the initial value ofα in is commensurate for α to start a sinusoidal form. Indeed, the initial values oḟ α at all the earlier visits to AB branch keep decreasing until it reaches a value that is consistent to begin the sinusoidal growth. Once this is satisfied, the monotonic increasing behavior from a to b is seen. As we will show this is the mechanism operating for high I and V case.
Case III, High I and V For this case, even on the AB branch,v/R cannot be ignored in Eq. (15) and thus one needs to use coupled Eqs. (21) and (15). Calculations follow much the same lines and give correct values for ∆t and ∆F on both the branches during the rapid jumps.
Again, we need to answer when exactly does the system know to switch from a rapid jumping mode to monotonically increasing on AB or decreasing mode on CD?
Consider the last of the rapid jumps from CD to AB (just prior to the point a) in Fig. 5(c). The corresponding point in the α plot [ Fig. 5(b)] is shown on an expanded scale in the inset. From this figure, it is clear thatα has a positive slope at k, though of small magnitude while at a, it has a value -9.7. The latter is close to the natural (negative) slope of α when it begins the descending branch of the sinusoidal form. On the other hand, the slope of α is positive at k and hence will not allow the growth to change over from a jumping mode to the sinusoidal growth form for α. One can note that the slopes at points of all the earlier visits to AB [see Fig. 5(b) inset] keep decreasing till the slope becomes negative required for the monotonically decreasing trend of α. This is exactly the same mechanism for I = 10 −2 and V = 1 also, for the low v branch, except that in this case, even the sign of the slope is incompatible for all the points prior to a in Fig. 5 c. The mechanism operating on CD (i.e., at the switching from jumping mode to monotonically decreasing nature of F ) is essentially the same but arguments are a little more involved and hence they are not presented. Now, we consider the causes leading to the maximum and minimum values taken by F being much more than permitted by f (v, V ). As this is dominant for V = 4, we illustrate this using Fig. 5(a) and 5(c). We first note that Eq. (10) constrains the dynamically changing values of F (t) and α(t) to the stationary values of f (v, V ). Clearly, this implies that F = f (v, V )/(1 + sinα). A rough estimate of F max can be obtained by F max ∼ f max /(1+α min ) with F min determined by α max . This relation can be easily verified by using the numerical values of α. small I, α is small, we should not find much difference between F max (F min ) and that of f .

VI. SUMMARY AND CONCLUSIONS
We first summarize the results before making some relevant remarks. We have carried out a study of the dynamics of an adhesive roller tape using a differentialalgebraic scheme used for singular set of differential equations. The algorithm produces stick-slip jumps across the two dissipative branches as a consequence of the inherent dynamics. Our extensive simulations show that the dynamics is much richer than anticipated earlier. In particular the influence of inertia is shown to be dramatic. For instance, even at low inertia, for small values of V , the influence of inertia manifests with jumps of the orbit occurring even before F reaches f max (or f min ) which is quite unexpected. More dominant is its influence for high I both for low V and high V , though it is striking for the latter case. Following the reasoning used in the PLC effect, we introduce a dynamized curve f (v, V ) as resulting from competing time scales of internal relaxation and imposed pull speed. The modified peel force function leads to the decreasing trend in the magnitude of ∆F with increasing pull velocity, a feature observed in experiments. We have also recaptured the essential features of the dynamics by a set of approximations valid in different regimes of the parameter space. These approximate solutions illustrate the influence of various time scales such as that due to inertia, the elasticity of the tape and that determined by the stationary peel force f (v, V ). We also find the unusual canard type of solutions.
Here, it is worthwhile to comment on the dynamical features of the model. The numerical results themselves are too complex to understand. A striking example of this is the series of force drops seen on the descending branch of the pull force [ Fig. 4(c)]. This result is hard to understand as it would amount to a partial relaxation of the pull force. However, a partial relaxation is only possible in the presence of another competing time scale (other than the imposed time scale). Another example is the jumping of the orbit for low I case, from AB to CB and vice versa even before the pull force reaches the extremum values of f (v, V ). For this reason, we have undertaken to make this complex dynamics transparent using a set of approximations. The basic idea here is to solve a single equation (or at most two equations as in the high I and V case) which incorporates all the relevant time scales. This method not only captures all the results to within 10% error but it also clearly brings out the regimes of parameter space where these time scales become important. This analysis also shows that the time scale due to inertia of the roller tape shows up even for low I which comes as a surprise as one expects that for low inertia, the orbit should stick to the stationary peel function. (Recall that for low inertia, equations have been approximated by Lienard type of equations by Maugis and Barquins [1].) Our approximate equations demonstrate that a crucial role in inducing the jumps even at low inertia is played by the high frequency oscillations resulting from the inertia of the roller tape. As for high inertia (both for low and high pull velocities), the time scale due to inertia is responsible for the partial relaxation of F as shown.
A few comments may be in order on the bursting type of oscillations in the peel velocity. Bursting type of oscillatory behavior are commonly seen in neuro-biological systems [25]. Conventionally, bursting type oscillations arise in the presence of homoclinic orbit [25]. Such bursting type of oscillations have also been modeled using one dimensional map [26]. However, it is clear that the mechanism for bursting type of oscillations in our case is different. In our case, this arises due to the fact that the orbit is forced to jump between the stable manifolds as a result of competing time scale of inertia and the time scale for the evolution of v. (We note that the latter itself includes more than one time scale [see Eq. (21)], namely the contribution from the slopes of the stable parts of the stationary curve f (v, V ) and that due to elasticity of the tape. ) The bunching of the spikes in v is the result of f (v, V ) becoming flat for large V and I. One other comment relates to canard type solutions. Figure 7 shows one such solution. As mentioned, these type solutions arise from sticking to the unstable manifold. In fact, a similar type of solution is seen in Fig. 5(c). As noted earlier, all the jumps from CD to AB or vice versa always occur when the peel velocity reaches the limiting value where f ′ (v, V ) = 0. However, it can be seen from this figure, the orbit starting from c monotonically decreases well into the unstable part of f . Thus, this solution also has the features of canards. It must be stated that our approximate solutions cannot capture the behavior of canards.
Finally, the results presented in this paper are on the nature of dynamics of the model equations which so far had defied solution. However, comparison with exper-iments has been minimal largely due to the paucity of quantitative experimental findings as stated earlier. Our analysis shows that the model predicts periodic, sawtooth [1], as well as chaotic solutions as reported in [3]. The high magnitude of the Lyapunov exponents for the chaotic solutions in the low pull velocities is consistent with that reported earlier [3]. We note that the other quantitative experimental feature reported by Refs. [1,3] is the decreasing trend of the average force drop magnitudes as a function of the pull velocity is also captured by our model (Fig. 6), a result that holds for both low and high inertia. This result is a direct consequence of the dynamization of the peel force function, i.e., dependence of the peel force on the pull velocity. We note here that the complex dynamics at high velocities (see Fig. 5) is a direct result of the unstable part of dynamized curve, f (v, V ), shrinking to zero. To the best of our knowledge, this is first time the result in Fig. 6 has been explained. As the hypothesis of dynamization captures the decreasing trend of the force drops, it also suggests that the underlying mechanism of competing time scales responsible for the peel force depending on the pull velocity is likely to be correct as in the PLC effect. Clearly, a rigorous derivation of the peel force function from microscopic considerations that includes the effect of the viscoelastic glue at the contact point is needed to understand the dynamics appropriately.