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.

## I. INTRODUCTION

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 *T*_{up}. Another string is attached at the bottom of the ball, and its tension is denoted as *T*_{down}. 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 obtain^{3,4}

where $ \omega 0 \u2261 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, *T _{c}*. According to Eq. (1), the lower string will break if

*α*>

*α*≡

_{c}*ω*

_{0}

*π*

^{−1}

*T*, 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 $ \alpha \u2208 [ 1 / 3 \alpha c , 1 / 2 \alpha c ] $. In this sense, which string breaks is not really a matter of pull or jerk.

_{c}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 capacity^{6} and thermal conductivity^{7–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 model^{15–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}

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.

## II. MODEL

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 $ \omega 0 \u2261 K / m $. The number of beads is *N*, and their equilibrium positions in the absence of external force is denoted by *x _{j}* =

*jl*, where

*l*means the equilibrium length of the spring (

*j*= 1,…,

*N*). The total length of the system in equilibrium therefore equals

*x*=

_{N}*Nl*≡

*L*. The displacement of the

*j*th bead from its equilibrium position is denoted by

*y*. The number of springs is also

_{j}*N*, and the extension of the

*j*th spring is

*z*≡

_{j}*y*−

_{j}*y*

_{j}_{−1}. The

*N*th 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:

where Γ is a friction coefficient. With the Kronecker delta *δ _{k}*

_{,}

_{l}and two auxiliary variables

*y*

_{0}≡ 0 and

*y*

_{N}_{+1}≡

*y*, all the above cases can be covered by the following expression:

_{N}where *j* = 1,…, *N*. In this notation, *y*_{0} = 0 and *y _{N}*

_{+1}=

*y*can be regarded as boundary conditions of Eq. (3). In particular,

_{N}*y*

_{0}= 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.,

*y*= 0 and

_{j}*dy*/

_{j}*dt*= 0 for every

*j*at

*t*= 0. In a dimensionless form, the dynamics is now rewritten as

where *τ* ≡ *ω*_{0}*t*, *ψ _{j}* ≡

*y*/

_{j}*l*,

*γ*≡ Γ/(2

*mω*

_{0}), 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

*ψ*is decomposed into homogeneous and particular parts, denoted by $ \psi j ( h ) $ and $ \psi j ( p ) $, respectively, to have $ \psi j = \psi j ( h ) + \u2009 \psi j ( p ) $. It is easy to see that $ \psi j ( p ) = j $ constitutes a particular solution for

_{j}*j*= 1,…,

*N*, with $ \psi 0 ( p ) \u2261 0 $ and $ \psi N + 1 ( p ) \u2261 \psi N ( p ) $. On the other hand, $ \psi j ( h ) $ satisfies the following homogeneous equation:

with $ \psi 0 ( h ) \u2261 0 $ and $ \psi N + 1 ( h ) \u2261 \psi N ( h ) $. The initial conditions are $ \psi j ( h ) = \u2212 j $ and $ d \psi j ( h ) / d \tau = 0 $ for *j* = 1,…, *N* at *τ* = 0. To construct a solution, one has to choose $ \psi j ( h ) \u221d \u2009 sin \u2009 k j $, considering $ \psi 0 ( h ) = 0 $. The other boundary condition $ \psi N + 1 ( h ) = \psi N ( h ) $ is then rewritten as $ sin \u2009 k ( N + 1 ) = sin \u2009 k N $, which quantizes the wavenumber as $ k n = ( n + 1 / 2 ) \pi / ( 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 ) \u2211 j = 1 N \u2009 sin \u2009 k n j \u2009 sin \u2009 k m j = \delta m n $. With these basis functions, the homogeneous solution is represented as $ \psi j ( h ) = \u2211 n = 0 N \u2212 1 a n ( \tau ) \u2009 sin \u2009 k n j $. Substituting this into Eq. (5), one sees that the coefficients have to satisfy

with $ \Omega n \u2261 2 \u2009 sin ( k n / 2 ) $, which is just the dispersion relation for phonons. The above differential equation can be solved by $ a n ( \tau ) = A e \mu n + \tau + B e \mu n \u2212 \tau $ with $ \mu n \xb1 \u2261 \u2212 \gamma \xb1 \gamma 2 \u2212 \Omega n 2 $ and arbitrary constants *A* and *B*. The constants are determined by applying the orthogonality relation to the initial conditions as follows:

where $ c n \u2261 \u2009 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:

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

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

It turns out that the case of *N* ≫ 1 greatly simplifies the analysis, making $ k n \u2248 \kappa n \u2261 ( n + 1 / 2 ) \pi / N , \u2009 \Omega n \u2248 \kappa n $, and $ \mu n \xb1 \u2248 \eta n \xb1 \u2261 \u2212 \gamma \xb1 \gamma 2 \u2212 \kappa n 2 $ for finite *n*. If a continuous variable *ξ* ≡ *x*/*l* is introduced to replace the integer index *j*, the Green's function becomes

with the displacement field,

as depicted in Fig. 3(a). The rescaled extension of the *j*th spring, $ \varphi j \u2261 \psi j \u2212 \psi j \u2212 1 = z j / l $, is approximated by $ \varphi ( \xi , \tau ) \u2261 ( \u2202 / \u2202 \xi ) \psi ( \xi , \tau ) $ evaluated at *ξ* = *j*. In terms of Eq. (13), it is written as

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

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

by calculating geometric series. When the argument *r* is away from (2*p* + 1)*N* for any integer *p*, the magnitude of *H _{N}*(

*r*) is small because of

*N*in the denominator. If

*r*− (2

*p*+ 1)

*N*=

*ϵ*≪ 1, on the other hand, the cosine in the denominator behaves linearly as

*ϵ*varies. It implies that

*H*is well approximated by the normalized sinc function, $ sinc \pi ( \u03f5 ) \u2261 \u2009 sin \u2009 \pi \u03f5 / ( \pi \u03f5 ) $, and the sign depends on

_{N}*p*as follows:

To sum up, *H _{N}*(

*r*) can be regarded as a train of sinc-typed impulses at

*r*= (2

*p*+ 1)

*N*. The factor of (−1)

^{p}means that two neighboring impulses have different signs, so the period of

*H*is 4

_{N}*N*in total. It would thus be useful to consider a convoluted function $ W ( r ) \u2261 sinc \pi ( r ) * \u0428 4 N ( r ) $, where $ \u0428 4 N ( r ) $ is the Dirac comb with periodicity of 4

*N*. The alternating impulse train is then described as $ H N ( r ) \u2248 W ( r \u2212 N ) \u2212 W ( r \u2212 3 N ) $ to a good approximation. Furthermore, one may simply take $ W ( r ) \u2248 \u0428 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 ) \u2248 \u0428 4 N ( r \u2212 N ) \u2212 \u0428 4 N ( r \u2212 3 N ) $. Now, the kernel function simplifies to

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

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 *j*th spring as follows:

where each of $ \nu max \xb1 $ 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 ( \tau ) = sin ( 2 \pi \tau / \tau 0 ) $ with *τ*_{0} = 4*N*. As expected, this induces resonant behavior [Fig. 3(c)] because *∂G*/*∂ξ* has 4*N*-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)].

## III. APPLICATION

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

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

*ξ*

At *ξ* = *N*, it reduces to

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 4*N* 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:

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 *τ* = 2*N*, 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

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

*ψ*and

_{j}*β*by a common factor

*λ*> 0

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

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 4

*N*. 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

*j*th spring is thus written as

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}

*ϕ*, it is a harmonic chain defined by

_{j}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

*N*th spring will break when

*β*is large because it is the earliest one that extends to

*ζ*. Conversely, small

_{c}*β*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

*τ*= 3

*N*−

*ξ*

^{*}(mod 4

*N*). In theory, therefore, it is possible to break every spring all at once by choosing a suitable value of

*β*so that every

*ϕ*reaches the threshold

_{j}*ϕ*at the same time, which may happen at

_{c}*τ*= 2

*N*. 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 *N*th 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

*ϕ*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

_{N}*γ*, the one that breaks will always be the

*N*th spring where the force is acting.

## IV. DISCUSSION AND CONCLUSION

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 breaking^{3–5} is still a theoretical possibility in this many-body system when driven by a ramp force *F* ∝ *t*. 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 $ \Delta t \u223c O ( \omega 0 \u2212 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 Δ

*t*∝

*m*

^{1∕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}

## ACKNOWLEDGMENTS

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).

### APPENDIX: SAMPLE CODE

Here, I present a Python code to simulate the dynamics of *N* = 30 beads with the RK4 method.^{19} Note that it includes *y*_{0} = 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 *y _{j}* 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