This work investigates the dynamics of a one-dimensional homogeneous harmonic chain on a horizontal table. One end is anchored to a wall; the other (free) end is pulled by an external force. A Green's function is derived to calculate the response to a generic pulling force. As an example, I assume that the magnitude of the pulling force increases with time at a uniform rate β. If the number of beads and springs used to model the chain is large, the extension of each spring takes a simple closed form, which is a piecewise-linear function of time. Under an additional assumption that a spring breaks when its extension exceeds a certain threshold, results show that for large β the spring breaks near the pulling end, whereas the breaking point can be located close to the wall by choosing small β. More precisely, the breaking point moves back and forth along the chain as β decreases, which has been called “anomalous” breaking in the context of the pull-or-jerk experiment. Although the experiment has been explained in terms of inertia, its meaning can be fully captured by discussing the competition between intrinsic and extrinsic time scales of forced oscillation.

The pull-or-jerk (or “inertia ball”) experiment is commonly used to demonstrate Newton's laws of motion: As depicted in Fig. 1, a ball of mass m is hung from a ceiling by a string, whose tension is denoted as Tup. Another string is attached at the bottom of the ball, and its tension is denoted as Tdown. A downward “jerky” force F is exerted at the end of the lower string, and the question is which string breaks. We know from our experience that the answer depends on how quickly the magnitude of the force changes: If the force suddenly increases, the lower string breaks. If, on the other hand, the force increases slowly, then the upper string breaks. This phenomenon is sometimes explained as due to the inertia of the ball, described as a “tendency to resist changes in motion” (see, e.g., Refs. 1 and 2). However, this does not clearly answer why the outcome varies with the jerkiness of the driving force, if all we know about mass is that it is a constant of motion independent of any specific dynamic process. Moreover, this experiment has a counter-intuitive aspect, which may be completely baffling unless one makes good sense of inertia: Suppose that the force has constant jerkiness α so that F = αt for t > 0. By approximating each string as a spring with force constant K, one can solve the equation of motion to obtain3,4

{ T up = α t α ω 0 1 sin ω 0 t , T down = α t ,
(1)

where ω 0 K / m and the gravitational force has been neglected. Interested readers are referred to Ref. 5 for a thorough analysis of this system. A simple assumption on the failure behavior is that the string breaks when T exceeds a certain threshold, say, Tc. According to Eq. (1), the lower string will break if α > αc ≡ ω0π−1Tc, which is consistent with our understanding. At the same time, it also predicts the existence of “anomalous” breaking, which means that the force can break the lower string with even smaller jerkiness α [ 1 / 3 α c , 1 / 2 α c ] . In this sense, which string breaks is not really a matter of pull or jerk.

Fig. 1.

Typical configuration of the pull-or-jerk experiment consisting of two elastic strings and a ball of mass m. The upper string has tension Tup. The lower one has Tdown and a time-dependent force F is exerted at its end.

Fig. 1.

Typical configuration of the pull-or-jerk experiment consisting of two elastic strings and a ball of mass m. The upper string has tension Tup. The lower one has Tdown and a time-dependent force F is exerted at its end.

Close modal

A natural question would be how this analysis generalizes to a chain of many beads and springs,5 which I wish to address in this work. A harmonic chain is a useful starting point to investigate properties of a macroscopic system near equilibrium. It has been used to understand basic statistical properties of solids such as heat capacity6 and thermal conductivity7–12 and the breaking strength of a polymer chain.13,14 A ladder of resistively and capacitively shunted Josephson junctions can also be approximated by a harmonic chain through the mapping to a locally coupled one-dimensional Kuramoto model15–17 when phase differences are small. Specifically, I will consider a harmonic chain consisting of identical beads and springs on a horizontal table as depicted in Fig. 2. One end of the chain is fixed to the wall, and the other end is driven by a time-varying external force F. The primary goal of this paper is to give students a precise picture of this general many-body pull-or-jerk experiment. Following Ref. 3, the spring is assumed to have a threshold of deformation above which it ceases to obey Hooke's law. The question is which spring is the earliest that reaches the threshold under a given external force, and this spring will be regarded as a breaking point of the chain. My finding is that the oscillatory motion in Eq. (1) manifests itself as a wave traveling across this many-body system, implying that one should consider two competing time scales, one for external driving and the other for internal wave dynamics, to understand the failure behavior. This study can also be thought of as an advanced exercise for physics majors because a harmonic chain is a representative mechanical example that is analytically soluble by means of undergraduate-level mathematics.18 

Fig. 2.

Schematic representation of a harmonic chain with N = 5 beads. Every bead has mass m, whereas the spring is massless. The spring constant is K for every spring. The leftmost spring is fastened to the wall whose displacement is fixed as y0 = 0, and the rightmost bead is driven by an external force F. The whole system is placed on a horizontal table so that the gravitational force does not enter the equation of motion.

Fig. 2.

Schematic representation of a harmonic chain with N = 5 beads. Every bead has mass m, whereas the spring is massless. The spring constant is K for every spring. The leftmost spring is fastened to the wall whose displacement is fixed as y0 = 0, and the rightmost bead is driven by an external force F. The whole system is placed on a horizontal table so that the gravitational force does not enter the equation of motion.

Close modal

This work is organized as follows: In Sec. II, I calculate the Green's function for the model system by solving the full equation of motion. It gives an approximate formula which holds in the continuum limit. The case of a linearly increasing force is then investigated in Sec. III under the assumption that friction is negligibly small. The analytic result is compared with numerical integration of the equations of motion. I discuss implications of the observed behavior and conclude this work in Sec. IV. A sample Python code is provided in the  Appendix.

Consider longitudinal waves on a harmonic chain consisting of beads and springs. Each bead has mass m and every spring has the same spring constant K, and the square root of their ratio is defined as ω 0 K / m . The number of beads is N, and their equilibrium positions in the absence of external force is denoted by xj = jl, where l means the equilibrium length of the spring (j = 1,…, N). The total length of the system in equilibrium therefore equals xN = Nl ≡ L. The displacement of the jth bead from its equilibrium position is denoted by yj. The number of springs is also N, and the extension of the jth spring is zj ≡ yj − yj−1. The Nth bead is pulled by a time-dependent external force F(t), where t denotes time. The equations of motion can thus be written as follows:

{ m d 2 d t 2 y 1 = K ( 2 y 1 + y 2 ) Γ d d t y 1 m d 2 d t 2 y j = K ( y j 1 2 y j + y j + 1 ) Γ d d t y j for 1 < j < N m d 2 d t 2 y N = K ( y N 1 y N ) Γ d d t y N + F ( t ) ,
(2)

where Γ is a friction coefficient. With the Kronecker delta δk,l and two auxiliary variables y0 ≡ 0 and yN+1 ≡ yN, all the above cases can be covered by the following expression:

m d 2 d t 2 y j = K ( y j 1 2 y j + y j + 1 ) Γ d d t y j + F ( t ) δ j , N ,
(3)

where j = 1,…, N. In this notation, y0 = 0 and yN+1 = yN can be regarded as boundary conditions of Eq. (3). In particular, y0 = 0 has direct physical meaning because the wall can be regarded as a fictitious bead with zero displacement (see Fig. 2). In addition, the system is initially at rest with zero displacements, i.e., yj = 0 and dyj/dt = 0 for every j at t = 0. In a dimensionless form, the dynamics is now rewritten as

d 2 d τ 2 ψ j = ( ψ j 1 2 ψ j + ψ j + 1 ) 2 γ d d τ ψ j + f ( τ ) δ j , N ,
(4)

where τ ≡ ω0t, ψj ≡ yj/l, γ ≡ Γ/(20), and f ≡ F/(Kl). It is convenient to choose f(τ) = u(τ) and solve Eq. (4) for τ > 0, where u means the Heaviside step function. The unknown ψj is decomposed into homogeneous and particular parts, denoted by ψ j ( h ) and ψ j ( p ) , respectively, to have ψ j = ψ j ( h ) + ψ j ( p ) . It is easy to see that ψ j ( p ) = j constitutes a particular solution for j = 1,…, N, with ψ 0 ( p ) 0 and ψ N + 1 ( p ) ψ N ( p ) . On the other hand, ψ j ( h ) satisfies the following homogeneous equation:

d 2 d τ 2 ψ j ( h ) = ψ j 1 ( h ) 2 ψ j ( h ) + ψ j + 1 ( h ) 2 γ d d τ ψ j ( h ) ,
(5)

with ψ 0 ( h ) 0 and ψ N + 1 ( h ) ψ N ( h ) . The initial conditions are ψ j ( h ) = j and d ψ j ( h ) / d τ = 0 for j = 1,…, N at τ = 0. To construct a solution, one has to choose ψ j ( h ) sin k j , considering ψ 0 ( h ) = 0 . The other boundary condition ψ N + 1 ( h ) = ψ N ( h ) is then rewritten as sin k ( N + 1 ) = sin k N , which quantizes the wavenumber as k n = ( n + 1 / 2 ) π / ( N + 1 / 2 ) , with n = 0, 1,…, N − 1. Note that the resulting basis functions are orthogonal in the sense that 4 / ( 2 N + 1 ) j = 1 N sin k n j sin k m j = δ m n . With these basis functions, the homogeneous solution is represented as ψ j ( h ) = n = 0 N 1 a n ( τ ) sin k n j . Substituting this into Eq. (5), one sees that the coefficients have to satisfy

d 2 d τ 2 a n + 2 γ d d τ a n = Ω n 2 a n ,
(6)

with Ω n 2 sin ( k n / 2 ) , which is just the dispersion relation for phonons. The above differential equation can be solved by a n ( τ ) = A e μ n + τ + B e μ n τ with μ n ± γ ± γ 2 Ω n 2 and arbitrary constants A and B. The constants are determined by applying the orthogonality relation to the initial conditions as follows:

a n ( 0 ) = A + B = 4 2 N + 1 j = 1 N ( j ) sin k n j = c n 2 Ω n 2 ,
(7)
d a n d τ ( 0 ) = A μ n + + B μ n = 4 2 N + 1 j = 1 N 0 × sin k n j = 0 ,
(8)

where c n sin [ ( 1 + N ) k n ] / ( N + 1 / 2 ) . After some algebra to compute A and B from the above set of equations, the solution is obtained as the following shifted discrete Fourier series:

ψ j ( h ) ( τ ) = n = 0 N 1 c n Ω n 2 γ 2 Ω n 2 ( μ n e μ n + τ μ n + e μ n τ ) sin k n j .
(9)

Using the connection between the Heaviside step function and the Dirac delta function, i.e., du/ = δ(τ), one readily obtains the Green's function as follows:

G j ( τ ) = d d τ ψ j ( τ ) = n = 0 N 1 c n γ 2 Ω n 2 ( e μ n + τ e μ n τ ) sin k n j ,
(10)

which is the response to f(τ) = δ(τ). Given any f(τ), the response can thus be calculated from the following convolution formula:

ψ j ( τ ) = 0 τ G j ( τ ) f ( τ τ ) d τ .
(11)

It turns out that the case of N ≫ 1 greatly simplifies the analysis, making k n κ n ( n + 1 / 2 ) π / N , Ω n κ n , and μ n ± η n ± γ ± γ 2 κ n 2 for finite n. If a continuous variable ξ ≡ x/l is introduced to replace the integer index j, the Green's function becomes

G ( ξ , τ ) n = 0 N 1 ( 1 ) n N γ 2 κ n 2 ( e η n + τ e η n τ ) sin κ n ξ ,
(12)

with the displacement field,

ψ ( ξ , τ ) = 0 τ G ( ξ , τ ) f ( τ τ ) d τ ,
(13)

as depicted in Fig. 3(a). The rescaled extension of the jth spring, ϕ j ψ j ψ j 1 = z j / l , is approximated by ϕ ( ξ , τ ) ( / ξ ) ψ ( ξ , τ ) evaluated at ξ = j. In terms of Eq. (13), it is written as

ϕ ( ξ , τ ) = 0 τ ξ G ( ξ , τ ) f ( τ τ ) d τ .
(14)
Fig. 3.

(a) Approximate displacement field ψ(ξ, τ) in Eq. (13) in response to the Heaviside step function, f(τ) = u(τ), when N = 30 and γ = 10−2. As time goes by, it converges to the particular solution, ψ(p) (ξ, τ) = ξ. (b) Propagator of Eq. (14), ∂G/∂ξ, when friction is ignored by setting γ = 0. (c) Extensions of springs under f ( τ ) = sin ( 2 π τ / τ 0 ) with τ0 = 4N. Equation (4) is integrated by the fourth-order Runge-Kutta method (RK4) with N = 30 and γ = 0, and the result agrees well with Eq. (23). (d) ϕ(ξ, τ) when f(τ) = δ(τ) + δ(τ − N).

Fig. 3.

(a) Approximate displacement field ψ(ξ, τ) in Eq. (13) in response to the Heaviside step function, f(τ) = u(τ), when N = 30 and γ = 10−2. As time goes by, it converges to the particular solution, ψ(p) (ξ, τ) = ξ. (b) Propagator of Eq. (14), ∂G/∂ξ, when friction is ignored by setting γ = 0. (c) Extensions of springs under f ( τ ) = sin ( 2 π τ / τ 0 ) with τ0 = 4N. Equation (4) is integrated by the fourth-order Runge-Kutta method (RK4) with N = 30 and γ = 0, and the result agrees well with Eq. (23). (d) ϕ(ξ, τ) when f(τ) = δ(τ) + δ(τ − N).

Close modal

If γ = 0 in Eq. (12), the kernel function ∂G/∂ξ takes the following form:

G ξ ( ξ , τ ) n = 0 N 1 ( 1 ) n N 2 sin κ n τ cos κ n ξ
(15)
= n = 0 N 1 ( 1 ) n N [ sin κ n ( τ + ξ ) + sin κ n ( τ ξ ) ]
(16)

inside the physical region, i.e., τ > 0 and 0 < ξ < N. Here, one can verify the following equality:

H N ( r ) n = 0 N 1 ( 1 ) n N sin κ n r = ( 1 ) N 1 sin π r 2 N cos π r 2 N
(17)

by calculating geometric series. When the argument r is away from (2p + 1)N for any integer p, the magnitude of HN(r) is small because of N in the denominator. If r − (2p + 1)N = ϵ ≪ 1, on the other hand, the cosine in the denominator behaves linearly as ϵ varies. It implies that HN is well approximated by the normalized sinc function, sinc π ( ϵ ) sin π ϵ / ( π ϵ ) , and the sign depends on p as follows:

H N ( 1 ) p sinc π ( ϵ ) .
(18)

To sum up, HN(r) can be regarded as a train of sinc-typed impulses at r = (2p + 1)N. The factor of (−1)p means that two neighboring impulses have different signs, so the period of HN is 4N in total. It would thus be useful to consider a convoluted function W ( r ) sinc π ( r ) * Ш 4 N ( r ) , where Ш 4 N ( r ) is the Dirac comb with periodicity of 4N. The alternating impulse train is then described as H N ( r ) W ( r N ) W ( r 3 N ) to a good approximation. Furthermore, one may simply take W ( r ) Ш 4 N ( r ) because the convoluted sinc function only modifies the peak shape without changing the essential physics. This leads to a particularly handy formula, H N ( r ) Ш 4 N ( r N ) Ш 4 N ( r 3 N ) . Now, the kernel function simplifies to

G ξ ( ξ , τ ) H N ( τ + ξ ) + H N ( τ ξ )
(19)
Ш 4 N ( τ + ξ N ) Ш 4 N ( τ + ξ 3 N ) + Ш 4 N ( τ ξ N ) Ш 4 N ( τ ξ 3 N ) ,
(20)

if τ > 0 and 0 < ξ < N. If the Dirac comb is written as an explicit sum of delta peaks, this can also be expressed as

G ξ ( ξ , τ ) δ ( τ + ξ N ) δ ( τ + ξ 3 N ) + δ ( τ ξ N ) δ ( τ ξ 3 N ) + δ ( τ + ξ 5 N ) δ ( τ + ξ 7 N ) + δ ( τ ξ 5 N ) δ ( τ ξ 7 N ) +
(21)
= ν = 0 ( 1 ) ν δ [ τ + ξ ( 2 ν + 1 ) N ] + ν = 0 ( 1 ) ν δ [ τ ξ ( 2 ν + 1 ) N ] .
(22)

On the (ξ, τ) plane, it is basically a pulse propagating back and forth between the two ends of the chain [Fig. 3(b)]. The pulse undergoes a phase shift of π every time it hits the free end on the right. Note that such soliton-like motion is due to the continuum approximation, in which the wave speed is given independent of the wave number when γ = 0. In a finite-sized system, the pulse will eventually disperse. Plugging Eq. (22) into Eq. (14), one approximately obtains the extension of the jth spring as follows:

ϕ j ( τ ) ν = 0 ν max + ( 1 ) ν f [ τ + j ( 2 ν + 1 ) N ] + ν = 0 ν max ( 1 ) ν f [ τ j ( 2 ν + 1 ) N ] ,
(23)

where each of ν max ± is defined as the greatest integer that makes positive the argument of every function in the summation. An example is the periodic driving force f ( τ ) = sin ( 2 π τ / τ 0 ) with τ0 = 4N. As expected, this induces resonant behavior [Fig. 3(c)] because ∂G/∂ξ has 4N-periodicity in time. It is also clear that one can observe constructive or destructive interference at a specific spring by sending pulses with an appropriate time interval [Fig. 3(d)].

If γ = 0 and f(τ) = βτ with a constant slope β, Eq. (13) yields

ψ ( ξ , τ ) = β τ ξ 16 β N 2 π 3 n = 0 ( 1 ) n ( 2 n + 1 ) 3 sin [ ( n + 1 2 ) π ξ N ] sin [ ( n + 1 2 ) π τ N ] .
(24)

A direct way to obtain Eq. (24) is to note that ψj = βτj forms a particular solution for Eq. (4) when γ = 0. For general γ > 0, however, one should employ the method of Green's functions. The displacement field is obtained by differentiating Eq. (24) with respect to ξ

ϕ ( ξ , τ ) = β τ 8 β N π 2 n = 0 ( 1 ) n ( 2 n + 1 ) 2 sin [ ( n + 1 2 ) π τ N ] cos [ ( n + 1 2 ) π ξ N ] .
(25)

At ξ = N, it reduces to

ϕ N ( τ ) ϕ ( ξ = N , τ ) = β τ ,
(26)

which means that the rightmost spring extends linearly in time. At the other end of the chain, i.e., at ξ = 0, the Fourier series on the right-hand side (RHS) of Eq. (25) describes a triangle wave of period 4N and amplitude βN, which behaves as −βτ between τ = 0 and N. Combining this result with the first term on the RHS of Eq. (25), one can see that the extension of the leftmost spring fastened to the wall is described by the following piecewise linear function:

ϕ 1 ( τ ) { 0 if 0 τ < N 2 β [ τ ( 2 ν + 1 ) N ] if ( 4 ν + 1 ) N τ < ( 4 ν + 3 ) N 4 β ( ν + 1 ) N if ( 4 ν + 3 ) N τ < ( 4 ν + 5 ) N ,
(27)

where ν = 0, 1, 2,…. It is plausible that the spring remains at rest when 0 ≤ τ < N because it takes time for the external perturbation to be transferred through N intermediate springs. However, the subsequent motion is not so self-evident: The spring suddenly begins to expand twice as fast as the rightmost one until the expansion stops abruptly at τ = 2N, and this pattern continues periodically. Application of Eq. (23) actually shows that every spring has such a discontinuity in the time derivative of ϕj except for j = N: From the shape of ∂G/∂ξ in Fig. 3(b), it is easily seen that

ϕ j ( τ ) { 0 if 0 τ < N j β ( τ + j N ) if ( 4 ν + 1 ) N j τ < ( 4 ν + 1 ) N + j 2 β [ τ ( 2 ν + 1 ) N ] if ( 4 ν + 1 ) N + j τ < ( 4 ν + 3 ) N j β ( τ j + N ) if ( 4 ν + 3 ) N j τ < ( 4 ν + 3 ) N + j 4 β ( ν + 1 ) N if ( 4 ν + 3 ) N + j τ < ( 4 ν + 5 ) N j ,
(28)

with ν = 0, 1, 2,…[see, e.g., ϕN∕2 in Fig. 4(a)]. Note that Eq. (28) is directly proportional to β. It is because Eq. (4) is linear and thus invariant under rescaling every ψj and β by a common factor λ > 0

λ d 2 d τ 2 ψ j = λ ( ψ j 1 2 ψ j + ψ j + 1 ) 2 γ λ d d τ ψ j + λ β τ δ j , N .
(29)
Fig. 4.

(a) Comparison between numerical results (points) and the approximate solutions in Eqs. (25)–(27) (lines), when γ = 0. The data points are obtained from the RK4 method with N = 30 and β = 1. Note that the extension of an intermediate spring such as ϕ15 lies between ϕ1 and ϕ30 all the time. (b) Rescaled mechanical energy of the chain [Eq. (30)], numerically obtained with the same set of parameters. (c) Which part breaks if every spring behaves as in Eqs. (31) and (32). Three different values of the threshold ζc are represented by the horizontal arrows. For large β, the corresponding ζc is small, so it is the Nth spring that extends to the threshold at the smallest τ (bottom arrow). For small β, the breaking point can be located closer to the wall because ζN is still below the threshold while others exceed it (middle arrow). However, the Nth spring can break with even smaller β (topmost arrow), which is called “anomalous” breaking. (d) Effects of friction when γ = 10−2 in Eq. (4). The other parameters are the same as in panel (a). Because of friction, the triangle wave of ϕ1 decays as time goes by, whereas ϕN is hardly affected.

Fig. 4.

(a) Comparison between numerical results (points) and the approximate solutions in Eqs. (25)–(27) (lines), when γ = 0. The data points are obtained from the RK4 method with N = 30 and β = 1. Note that the extension of an intermediate spring such as ϕ15 lies between ϕ1 and ϕ30 all the time. (b) Rescaled mechanical energy of the chain [Eq. (30)], numerically obtained with the same set of parameters. (c) Which part breaks if every spring behaves as in Eqs. (31) and (32). Three different values of the threshold ζc are represented by the horizontal arrows. For large β, the corresponding ζc is small, so it is the Nth spring that extends to the threshold at the smallest τ (bottom arrow). For small β, the breaking point can be located closer to the wall because ζN is still below the threshold while others exceed it (middle arrow). However, the Nth spring can break with even smaller β (topmost arrow), which is called “anomalous” breaking. (d) Effects of friction when γ = 10−2 in Eq. (4). The other parameters are the same as in panel (a). Because of friction, the triangle wave of ϕ1 decays as time goes by, whereas ϕN is hardly affected.

Close modal

In words, the sole effect of choosing a different value for β is to change the overall length scale. No matter how slowly the end of the harmonic chain is pulled, the periodic discontinuity will not disappear. Note also that Eqs. (26) and (27) provide envelopes for every ϕj in between [Eq. (28)], although the lines can sometimes coincide. As demonstrated in Fig. 4(a), the analytic predictions of Eq. (28) are well substantiated by direct numerical integration of Eq. (4). One may also check the mechanical energy per particle

ε = 1 N j = 1 N 1 2 ( 1 N d ψ j d τ ) 2 + 1 N j = 1 N 1 2 ( ψ j ψ j 1 N ) 2 ,
(30)

where the first and second terms represent the kinetic and potential parts, respectively [Fig. 4(b)]. The factor of 1/N inside each pair of parentheses is due to the fact that ψj ∼ O(N) [see Fig. 4(a)]. The kinetic part of Eq. (30) turns out to be a periodic function of τ with a period of 4N. The potential energy, on the other hand, keeps increasing as a quadratic function of τ.

Recall the assumption that every spring obeys Hooke's law up to some threshold ϕc > 0, above which the spring breaks.3–5 The restoring force of the jth spring is thus written as

f j res = { K ϕ j if | ϕ j | < ϕ c 0 otherwise ( i . e . , broken ) .
(31)

Then, the above calculation implies that the magnitude of β is an important factor to determine which spring breaks. It is related to the fact that a different value of β just rescales every ϕj with exactly the same factor. If λ = β−1 in Eq. (29) and ζj ≡ β−1ϕj, it is a harmonic chain defined by

d 2 d τ 2 ζ j = ( ζ j 1 2 ζ j + ζ j + 1 ) 2 γ d d τ ζ j + τ δ j , N ,
(32)

in which the threshold of a spring becomes ζc ≡ ϕc/β. Large jerkiness therefore maps to a low threshold in this derived system. As illustrated in Fig. 4(c), the Nth spring will break when β is large because it is the earliest one that extends to ζc. Conversely, small β can break a spring close to the wall. Precisely speaking, one can only specify the range of springs to break in the latter case. According to Eq. (28), all the springs between ξ = 0 and ξ*(≤N) are the most extended ones in this chain when τ = 3N − ξ* (mod 4N). In theory, therefore, it is possible to break every spring all at once by choosing a suitable value of β so that every ϕj reaches the threshold ϕc at the same time, which may happen at τ = 2N. In practice, however, this would mean that the breaking point becomes very sensitive to experimental noise and mechanical defects.

Another point of Fig. 4(c) is that one can break the Nth spring by pulling the end even more slowly, which proves the existence of “anomalous” breaking in this system.3–5 If this anomaly is hardly observed, the reason could be that friction is not negligible in any experimental situation.3 If γ is positive yet so small that only the second summation contributes in Eq. (12), the triangle wave in ϕ1 will gradually decay as indicated by eγτ in the summand. On the other hand, it is reasonable to guess that ϕN will not experience any notable change, considering that ∂G/∂ξ containing the friction term identically vanishes at ξ = N due to the boundary condition. Consequently, ϕ1 is expected to lie below ϕN in the long run [Fig. 4(d)]. It implies that the anomaly can indeed be diminished by friction, but the price is that it also becomes hard to locate the breaking point close to the wall. For sufficiently large γ, the one that breaks will always be the Nth spring where the force is acting.

To summarize, I have investigated the dynamics of a harmonic chain which is anchored to a wall at one end and subject to an external force at the other end. By using the method of Green's functions, one can calculate the response of the system to a general time-varying force. A simple expression is obtained when the system becomes a continuous medium composed of a large number of beads and springs [Eq. (23)]. With the simple failure behavior assumed in Eq. (31), anomalous breaking3–5 is still a theoretical possibility in this many-body system when driven by a ramp force Ft. A nontrivial difference from the common pull-or-jerk experiment is that every spring except the last one exhibits distinct stop-and-go behavior in its extension ϕj [compare Eq. (1) and Eq. (28)]. It implies that it roughly takes Δ t O ( ω 0 1 ) for the external perturbation to travel across a spring, and this is a fast process compared to system-wide dynamics when N is large. When one talks about the pull-or-jerk experiment in the context of Newton's law of inertia, the precise meaning is that Δtm1∕2. If the time is not enough to send an amount of energy across the chain, therefore, it will be the rightmost spring that breaks. In other words, there are two competing time scales: One is the intrinsic time scale of the chain, and the other is that of the driving force. The point is that the experiment should be understood in terms of these time scales of forced oscillation, in addition to the law of inertia.

From a technical point of view, the chain is described by a set of coupled, linear, ordinary differential equations. Although it looks much more difficult than the one-body counterpart as in Fig. 1, the problem can readily be handled by standard techniques such as separation of variables and Green's functions.18 It is instructive to check the validity of the analytic solution by performing numerical simulations, e.g., with the RK4 method as we have done throughout this work. Figure 4(a) has already shown consistency between the analytical and numerical approaches, but the agreement is actually striking in every detail, as demonstrated in Fig. 5. For reference, a sample Python code is provided in the  Appendix.19 Readers are also encouraged to extend the model to inhomogeneous or anharmonic cases to incorporate more realistic aspects of the chain dynamics.20 

Fig. 5.

Comparison between the analytical solution [Eq. (11)] and the numerical result obtained by the RK4 method when N = 30, γ = 0, and f(τ) = τ. Inset: Zoomed-in view of the small oscillatory part. The numerical data points are exactly on top of the analytical solution.

Fig. 5.

Comparison between the analytical solution [Eq. (11)] and the numerical result obtained by the RK4 method when N = 30, γ = 0, and f(τ) = τ. Inset: Zoomed-in view of the small oscillatory part. The numerical data points are exactly on top of the analytical solution.

Close modal

S.K.B. gratefully acknowledges discussions with Julian Lee, Hang-Hyun Jo, and Hiroyuki Shima. This work was supported by a research grant of Pukyong National University (2016).

Here, I present a Python code to simulate the dynamics of N = 30 beads with the RK4 method.19 Note that it includes y0 = 0 explicitly because the index of an array begins from zero by default. The array r contains both the displacement and velocity of every bead in such a way that its elements r[2*j] and r[2*j+1] correspond to yj and ( d / d t ) y j , respectively.

from _ _future_ _ import print_function, division # for Python 2

from numpy import empty, array, zeros

def drive(t):

return beta * t

def f(r,t):

 rdot = empty(2 * N1, float)

for i in range(0, 2 * N1, 2): # time derivative of displacement

  rdot[i] = r[i + 1]

for i in range(1, 2 * N1, 2): # time derivative of velocity

  if i ==1:

     accel = 0.

  elif i > 1 and i < 2 * N1−1:

               accel = r[i − 3] − 2 * r[i − 1]+ r[i + 1]−gamma * r[i]

  elif i ==2 * N1 − 1:

     accel = r[i − 3] − r[i − 1] − gamma * r[i] + drive(t)

  rdot[i] = accel

return rdot

start, end = 0., 150.       # time domain

max_step = 15000         # number of time steps

h = (end − start) / max_step # time increment

N1 = 31     # number of beads including the zeroth one (j = 0)

beta = 1.    # increasing rate of the driving force

gamma = 0. # friction coefficient

r = zeros(2 * N1, float)

for step in range(max_step):

 t = h * step

 k1 = h * f(r,t)

 k2 = h * f(r + 0.5 * k1,t + 0.5 * h)

 k3 = h * f(r + 0.5 * k2,t + 0.5 * h)

 k4 = h * f(r + 0.5 * k3,t + h)

 r + = (k1 + 2 * k2 + 2 * k3 + k4)/6

1.
P. G.
Hewitt
,
Conceptual Physics
, 11th ed. (
Addison-Wesley Longman
,
Reading, MA
,
2009
).
2.
Information is available in the “More inertia experiments” section of the Institute of Physics Web site, <http://practicalphysics.org/more-inertia-experiments.html>.
3.
M. A.
Heald
and
G. M.
Caplan
, “
Which string breaks?
,”
Phys. Teach.
34
(
8
),
504
507
(
1996
).
4.
G. M.
Caplan
and
M. A.
Heald
, “
Ye olde inertia demonstration
,”
Am, J. Phys.
72
(
7
),
860
862
(
2004
).
5.
H.
Shima
, “
Analytic expression for the pull-or-jerk experiment
,”
Eur. J. Phys.
35
(
6
),
065016
(
2014
).
6.
R. H.
Swendsen
,
An Introduction to Statistical Mechanics and Thermodynamics
(
Oxford U.P.
,
Oxford
,
2012
).
7.
Z.
Rieder
,
J. L.
Lebowitz
, and
E.
Lieb
, “
Properties of a harmonic crystal in a stationary nonequilibrium state
,”
J. Math. Phys.
8
(
5
),
1073
1078
(
1967
).
8.
J.
Florencio
, Jr.
and
M. H.
Lee
, “
Exact time evolution of a classical harmonic-oscillator chain
,”
Phys. Rev. A
31
(
5
),
3231
3236
(
1985
).
9.
M. H.
Lee
,
J.
Florencio
, Jr.
, and
J.
Hong
, “
Dynamic equivalence of a two-dimensional quantum electron gas and a classical harmonic oscillator chain with an impurity mass
,”
J. Phys. A: Math. Gen.
22
(
8
),
L331
L335
(
1989
).
10.
B.
Hu
,
B.
Li
, and
H.
Zhao
, “
Heat conduction in one-dimensional nonintegrable systems
,”
Phys. Rev. E
61
(
4
),
3828
3831
(
2000
).
11.
A.
Dhar
, “
Heat conduction in the disordered harmonic chain revisited
,”
Phys. Rev. Lett.
86
(
26
),
5882
5885
(
2001
).
12.
F.
Bonetto
,
J. L.
Lebowitz
, and
J.
Lukkarinen
, “
Fourier's law for a harmonic crystal with self-consistent stochastic reservoirs
,”
J. Stat. Phys.
116
(
1
),
783
813
(
2004
).
13.
T.
Doerr
and
P.
Taylor
, “
Breaking in polymer chains. I. The harmonic chain
,”
J. Chem. Phys.
101
(
11
),
10107
10117
(
1994
).
14.
C. F.
Lee
, “
Thermal breakage of a discrete one-dimensional string
,”
Phys. Rev. E
80
(3
),
031134
(
2009
).
15.
K.
Wiesenfeld
,
P.
Colet
, and
S. H.
Strogatz
, “
Synchronization transitions in a disordered Josephson series array
,”
Phys. Rev. Lett.
76
(
3
),
404
407
(
1996
).
16.
K.
Wiesenfeld
,
P.
Colet
, and
S. H.
Strogatz
, “
Frequency locking in Josephson arrays: Connection with the Kuramoto model
,”
Phys. Rev. E
57
(
2
),
1563
1569
(
1998
).
17.
B.
Daniels
,
S.
Dissanayake
, and
B.
Trees
, “
Synchronization of coupled rotators: Josephson junction ladders and the locally coupled Kuramoto model
,”
Phys. Rev. E
67
(
2
),
026216
(
2003
).
18.
M. L.
Boas
,
Mathematical Methods in the Physical Sciences
, 3rd ed. (
Wiley
,
Hoboken, NJ
,
2006
).
19.
M. E. J.
Newman
,
Computational Physics
(
CreateSpace Independent
,
United States
,
2013
).
20.
J.
Paturej
,
A.
Milchev
,
V. G.
Rostiashvili
, and
T. A.
Vilgis
, “
Polymer chain scission at constant tension-an example of force-induced collective behaviour
,”
EPL
94
(
4
),
48003
(
2011
).