Synchrony in networks of Franklin bells

The Franklin bell is an electro-mechanical oscillator that can generate a repeating chime in the presence of an electric field. Benjamin Franklin famously used it as a lightning detector. The chime arises from the impact of a metal ball on a metal bell. Thus, a network of Franklin bells can be regarded as a network of impact oscillators. Although the number of techniques for analysing impacting systems has grown in recent years, this has typically focused on low-dimensional systems and relatively little attention has been paid to networks. Here we redress this balance with a focus on synchronous oscillatory network states. We first study a single Franklin bell, showing how to construct periodic orbits and how to determine their linear stability and bifurcation. To cope with the non-smooth nature of the impacts we use saltation operators to develop the correct Floquet theory. We further introduce a new smoothing technique that circumvents the need for saltation and that recovers the saltation operators in some appropriate limit. We then consider the dynamics of a network of Franklin bells, showing how the master stability function approach can be adapted to treat the linear stability of the synchronous state for arbitrary network topologies. We use this to determine conditions for network induced instabilities. Direct numerical simulations are shown to be in excellent agreement with theoretical results.


Introduction
The history of the Franklin bell is long and well established. Although named after the American scientist Benjamin Franklin it was in fact invented by the Scottish Benedictine monk Andrew Gordon in Erfurt, Germany, around 1742. The bell converts electrical energy into mechanical energy in the form of a repeating mechanical motion and forms the basis for many modern day bell-chimes, from security alarms to school bells. Franklin made use of Gordon's idea by connecting one bell to his pointed lightning rod, attached to a chimney, and a second bell to the ground. One of his letters contains the following description (Franklin, 1962) In September 1752 I erected an Iron Rod to draw the Lightning down into my House, in order to make some Experiments on it, with two bells to give Notice when the Rod should be electrified. A contrivance obvious to every Electrician.
while the other bell conducts the discharge to the ground. The shuttling behaviour of the metallic ball between bells produces the chime and allows the detection of lightning. This is a prototypical example of an impact oscillator.
In recent years, a considerable amount of research has been devoted to the development of efficient techniques to analyse the dynamical behaviours of impacting systems. This has been motivated in part by challenges arising in control theory, population dynamics, chemistry, physics, biotechnologies, economics, industrial robotics, to name but a few (Samoilenko & Perestyuk, 1995;Yang, 2001;Catllá et al., 2008;di Bernardo et al., 2008a,b;Simpson & Kuske, 2018). Indeed, many real world systems can be characterized by instantaneous jumps or switches in behaviour, which may be created by impulsive interactions. In contrast to smooth dynamical systems, the analysis of such non-smooth systems is relatively underdeveloped. This is even more true at the network level. Thus, it is of interest to either adapt techniques from the theory of smooth dynamical systems or to develop entirely new ones. We do so here with a focus on synchronous periodic states in Franklin bell networks of arbitrary topology. For some demonstrations of Franklin bell networks we refer the reader to the growing number of videos that are being increasingly used in scientific outreach activities (RimstarOrg, 2012).
Synchrony is a common behaviour seen in networks of oscillators with graph Laplacian coupling (of which diffusive coupling is a classic example), and arises in many different areas of biology, engineering, ecology, and social sciences (Pikovsky et al., 2001;Nijmeijer & Rodriguez-Angeles, 2003;Pecora et al., 1997;Wang & Chen, 2002;Ariaratnam & Strogatz, 2001;Sorrentino et al., 2016;Pogromsky et al., 2002;Yu et al., 2014;Zhang et al., 2007;Arenas et al., 2008). In order to deal with the synchronization problem for coupled oscillators, Pecora and Carroll introduced the master stability function (MSF) technique for smooth systems (Pecora et al., 1997;Pecora & Carroll, 1998). Under some assumptions (identical oscillators and graph Laplacian coupling) they developed a network Floquet theory that can be diagonalized in the basis of the eigenvectors of the network connectivity matrix. This means that the stability of the synchronous orbit can be assessed in terms of a set of lower dimensional Floquet problems parameterized by the (possibly complex) eigenvalues of the network connectivity matrix. Recently, this method has been extended to treat diffusively coupled networks of non-smooth Filippov type (Filippov, 1988) and integrate-and-fire piecewise linear (PWL) oscillator models (Coombes & Thul, 2016;Nicks et al., 2018;Lai et al., 2018), making use of saltation operators. These have been widely used in the non-smooth dynamical systems community to treat the linearised evolution of small perturbations through switching manifolds (Müller, 1995).
As well as developing the mathematical techniques for handling a truly non-smooth Franklin bell network, we further introduce a new form of smoothing that circumvents the need for constructing saltation operators. At heart, this technique introduces a virtual linear dynamical system that smoothly connects the orbits before and after impact. The duration of this virtual trajectory (that bridges the impact) is chosen as a control parameter δ. In the limit that δ tends to zero the propagator for this virtual system recovers the saltation rule. Thus working with small but finite δ we may treat the Franklin bell network solely as a smooth system. As expected the MSF for the non-smooth network and the smoothed network show excellent agreement for small δ.
The organization of the paper is as follows. In Section 2, we provide a detailed description of a model for a single Franklin bell. We show how to construct periodic orbits and adapt Floquet theory using saltation operators to assess solution stability. We use this to determine the bifurcation diagram as a function of the restitution of the ball velocity upon impact with the bell. In Section 3, we present the new smoothing technique and show that it recovers the saltation operators previously constructed. Then in Section 4 we use the MSF technique to determine the stability of synchronous network states, for both smooth and non-smooth networks. Numerical examples are presented in Section 5 and shown Here a battery generates a constant voltage V across the circuit. An initially stationary ball, hanging midway between the two charged plates, will have a charge distribution that is positive on its right-hand side and negative on its left hand side. Left: an initial push to the left will cause a stationary ball to be attracted to the left-hand plate. Upon impact it will exchange charge with the plate and develop a net positive charge distribution. The ball and plate will then both have positive charge and repel each other. Right: the repulsive force from the left will cause an impact of the ball with the right-hand plate, where it can collect negative charge. The negatively charged metal ball will then be repelled from the negatively charged plate, and will move to the left. Thus, a repetitive impacting oscillation can develop. This is the basis for chiming in a Franklin bell. to be in excellent agreement with theory. Finally, in Section 6 we discuss the results in this paper and natural extensions of the work presented.

Model description and periodic orbits
In it simplest form the Franklin bell can be regarded as an electro-mechanical system consisting of two oppositely charged parallel plates (representing the metal bells) with a conductive particle (metal ball) which travels between them. This is suspended from an insulting wire hanging midway between the two plates. The polarization of the plates is maintained by a battery, such that the constant electric field between them generates an electrostatic force that causes the ball to move. Upon impact with a plate the ball reverses direction and moves toward the opposite plate. In this way a periodic impacting rhythm can be generated. This is illustrated in Fig. 1. To formulate a mathematical description of this idealized process we consider the metallic plates to be placed at the positions u = ±a, with u ∈ R. The ball that travels between the plates is governed by the dynamics of a forced damped simple harmonic (pendulum) oscillator. The direction of the forcing is determined by the sign of the charge carried by the ball at the instant before impact and is reversed after impact. The magnitude of this force is determined by the sum of the repelling and attracting electrostatic forces, and will be assumed to be a constant denoted by f . Thus, we are led to the equations of motion for a single Franklin bell as Here u,u, andü are the position, velocity and acceleration of the particle at time t, respectively. The damping coefficient is given by γ 1 > 0, γ 2 > 0 sets the natural frequency of the pendulum, and k ∈ R + is the coefficient of restitution upon impact. The impact times t i are determined implicitly from the conditions u(t i ) = ±a, i ∈ Z. The dramatic change in velocity at impact is governed by equation (2.2), whereu(t − i ) represents the velocity of the ball immediately before t = t i andu(t + i ) immediately thereafter, i.e.u(t ± i ) = lim t→t ± iu (t). It is clear from (2.1) that the right-hand side of the system changes discontinuously upon impact, so that it may be regarded as a Filippov system. Moreover, the system is impulsive because of the velocity jump at the impact times. Thus, we consider the basic Franklin bell model as a state-dependent impulsive system with discontinuous changes in the vector field at impact times. Models of this type exist in many in real world scenarios (Yang, 2001;Liu & Wang, 2006;Müller, 1995;Fredriksson & Nordmark, 2000;di Bernardo et al., 2001), and are exemplified by impact oscillators. Thus, it is natural to analyse the Franklin bell as a state-dependent impacting system (Bishop, 1994;Fredriksson & Nordmark, 2000).
It is first sensible to examine the fixed point structure of the model. Introducing v =u and denoting x = (u, v) T , then (2.1) and (2.2) can be written in a state-space form aṡ and g(x) = (u, −kv) T . As well as the jump rule g for describing what happens at impact it is convenient to introduce two indicator functions h = h ± (x) = u ± a that determine the times of impact according to h(x(t i )) = 0. From (2.3) the stability of any equilibrium points is determined by the eigenvalues of the matrix A. These are easily calculated as λ ± = (−γ 1 ± √ (γ 2 1 − 4γ 2 ))/2. Remembering that γ 1 and γ 2 are both positive, we see that if a fixed point exists then it is stable (being a node for γ 2 1 − 4γ 2 0 and a focus otherwise). Formally, equilibrium points can be calculated as (u, v) = (±f /γ 2 , 0). Consequently if |f | < γ 2 |a| then both fixed points will be between the two plates, and otherwise they will be virtual (lying outside of the physically accessible region). This latter case will guarantee the existence of impacts, and is the one we focus on for the rest of the paper since it is a necessary condition for the existence of periodic orbits, and hence chiming in a Franklin bell.

Construction and stability of periodic solution
In general it is very hard to find closed form solutions for periodic orbits in nonlinear dynamical systems. However, since (2.3) is a PWL system, it can be solved exactly in regions of phase space where v > 0 and v < 0, respectively, and solutions glued together to construct periodic orbits. Consider now a periodic motion that starts at t = 0 + at Plate 1 (see Fig. 1) and returns to the same point after a period Δ. Let us denote the time of flight for the trajectory from u = −a to u = a with v > 0 by Δ 1 . An explicit form for this trajectory can be constructed using matrix exponentials and initial data The impact time t 1 is determined by the condition u(t − 1 ) = a, and the time of flight is simply Δ 1 = t 1 . An application of the jump rule can then be used to determine new initial data for the trajectory in the region, where v < 0. This yields . (2.7) Denoting the time of flight for the trajectory from u = a to u = −a with v < 0 by Δ 2 , then the corresponding trajectory is The impact time t 2 is determined by the condition u(t − 2 ) = −a, and the time of flight is simply Δ 2 = t 2 − t 1 . An application of the jump rule at time t 2 then gives ( , and for the orbit to be periodic this must match the initial data (u(0 + ), v(0 + )) = (−a, v 0 ). Thus, a periodic orbit, parametrized by (t 1 , t 2 , v 0 ), will exist if there is a solution to the three simultaneous nonlinear algebraic equations If a solution exists then the period of oscillation for a periodic orbit with x(t) = x(t + Δ) is given by To determine the stability of periodic solutions we should be careful about the evolution of perturbations through the impacting manifolds (where u = ±a). If we denote a perturbation to the periodic orbit by δx(t), then between impacts the linearised evolution of small perturbations is governed by with solution δx(t) = exp(At)δx(0). The application of a saltation operator can then be used to map a perturbed trajectory across the impact manifold. The methodology for the construction of the relevant saltation operation is now well developed, see e.g. Leine & Nijmeijer (2004). In our case this saltation operation can be expressed in terms of a matrix and Dg is the Jacobian of g. Using the above we may construct the saltation matrices at t = t 1 (where u = +a) and t = t 2 (where u = −a) as K(t 1 ) and K(t 2 ) respectively, where (2.12) Thus, after one period of oscillation a perturbed trajectory will have evolved according to the formula Thus, the periodic orbit will be stable if the eigenvalues of Q lie within the unit disc. Since for a planar system one of the Floquet multipliers is equal to one (corresponding to tangential perturbations along the orbit) there is only one non-trivial eigenvalue of Q to consider. If we denote this by e σ Δ and use the result that det Q = e σ Δ × 1, then we have that (2.14) A periodic orbit will be stable if σ < 0. Thus, if k < 1 then all periodic orbits must be stable. However, if the coefficient of restitution were taken to be greater than one (corresponding to injecting energy into  For k 1 there is only one stable periodic orbit. At k = 1 an unstable orbit with infinitely large amplitude is born whose amplitude decreases with a further increase in k. For k > 1 stable and unstable periodic orbits co-exist until k 1.663 where they are lost in a saddle-node bifurcation of periodic orbits. The amplitude of the stable periodic orbit is always less than that of the unstable orbit. Black solid (dashed) and red solid (dashed) curves represent Floquet exponents and the amplitude of stable (unstable) periodic orbits, respectively. the system at impact) then it would be possible for unstable periodic orbits to exist. An example of a co-existing stable and unstable periodic orbit for k > 1 is shown in Fig. 2. Here a stable periodic orbit is encircled by a an unstable periodic orbit with a large amplitude.
A bifurcation diagram, summarizing the properties of periodic orbits under variation in the coefficient of restitution k is shown in Fig. 3. Here we see that for k < 1 there is only one stable periodic orbit, whilst for k > 1 a new unstable period orbit of large amplitude can be created. With increasing k, it is ultimately destroyed in a saddle-node bifurcation of periodic orbits. Mechanically, the case with k > 1 corresponds to energy being pumped into the system at impact, as in many pinball machines, and is often referred to as active impact (Pring & Budd, 2011).

A piecewise linear smoothing technique
Although the non-smooth system can be treated rigorously with the use of saltation operators it is of interest to consider a smoothed version of the model, which can be analysed with more traditional techniques. The choice of smoothing is somewhat arbitrary and one may consider a variety of approaches-a discussion can be found in Jeffrey (2018). If the model is written using potentials, then the non-smooth system has an infinitely steep potential at the two plates, which could be replaced by a potential function with finite but steep gradient at the plates. Instead here we choose to append new dynamical rules at the end plates (and the regions beyond them), remove the strict impact condition, and allow trajectories to cross through the switching manifolds. We now envisage trajectories, beyond the plates, that can smoothly connect to those within the plates. If the latter are determined by the original non-smooth system then this effectively provides a smoothing. Although preserving the shape of an orbit this does not preserve its proper duration as further time is needed to traverse the region outside the plates. If the time-of-flight could be reduced to zero outside the plates then this would recover the truly nonsmooth trajectory. Here we show that this can be achieved with a simple choice of linear dynamical system outside the plates.
The formal description of the smoothed model is written by augmenting the original model, given by (2.1) and (2.2), in the following way: for as yet undetermined matrices A R,L ∈ R 2×2 and vectors f R,L ∈ R 2 . Each of the two new linear dynamical system is thus determined by six unknown parameters (four for A R,L and two for f R,L ). These can be computed from matching conditions at the points, where u = ±a such that the solution for x andẋ is continuous and respects the rule for restitution. For example, if we denote the value of x when u = +a by This gives a total of six equations for six unknowns, parameterized by the time-of-flight δt. The equation for x(t) can be determined explicitly using matrix exponentials, as in (2.6) and (2.8), using the matrix A or A R,L as appropriate. The six simultaneous nonlinear equations can be solved numerically. Similarly, we may match at u = −a and obtain a similar system of equations (under the interchange of labels R to L). An illustration of this process is given in Fig. 4.
In Fig. 5 we show two examples of trajectories constructed by patching together matrix exponential solutions in the regions u < −a, |u| a, and u > a, subject to the smoothing process described above. As we take δt smaller and smaller we find that the smooth trajectory approaches that of the non-smooth where Δ 1 , Δ 3 are the times of flight in the region |u| a and Δ 2 = Δ 4 = δt in the regions where |u| > a. Moreover, the propagators in the regions u < −a and u > a, exp(A L t) and exp(A R t), respectively, approximate the relevant saltation matrices. The numerical evidence for this is provided in Fig. 6, where we compare the components of K(t 1 ) (see (2.12)) with those of exp(A R δt).
Although δt is under our control it is not guaranteed that this time-of-flight can be made arbitrarily small. Here, we have only provided numerical evidence that this is the case, and have not provided a formal proof. Rather we have presented a practical method for smoothing systems with hard impacts, obviating the need for the construction of saltation matrices. In the next section we show how to treat networks of interacting Franklin bells, with both hard and smoothed impacts.

A Franklin bell network
A Franklin bell network can easily be constructed by serial extension of the network shown in Fig. 1. One simply hangs more metal balls from a cross-bar and inserts a metal bell between each suspended ball. Other topologies are, of course, possible that leads us to the consideration of general Franklin bell networks. The vertices of such a network can be represented by the bell-ball-bell combination and network edges by the interactions between them. From a modelling perspective the interactions between nodes are mediated by the vibrations communicated through the cross-bar. This is very reminiscent of a system of Huygens clocks (Huygens, 1893), albeit where the clocks are impact oscillators rather than smooth limit-cycle oscillators. There is now a vast literature on the study of the latter, see e.g. Dörfler & Bullo (2014), though far less work has been done on networks of impact oscillators. The exception to this perhaps being the recent work of Shiroky & Gendelman (2016) who examined a linear array of Franklin bell oscillators. They analysed the properties of localized states (breathing modes), whereby only one of the network nodes made repetitive impacts. The stability and bifurcation of these localized states was determined using a Fourier-based Floquet theory adapted to cope with local impulsive (Dirac-delta) effects. Thorin et al. (2017) have also considered a similar problem, and most recently James et al. (2018) have highlighted some of the many open problems in the study of impact oscillator networks. Here we focus on globally periodic synchronous impacting behaviour and show how to augment techniques from the network science of smoothly coupled limit cycles to treat impact oscillators. Moreover, by exploiting the PWL nature of a Franklin bell network between impacts we show how to readily construct the MSF. This is a powerful tool for determining the stability of a synchronous orbit in a network of arbitrary topology.
We begin by describing the construction of the MSF for an impulsive nonsmooth Franklin bell network, and then indicate how to perform the same calculation for a smoothed system.

Master stability function for a non-smooth Franklin bell network
Consider an impacting Franklin bell network with N identical nodes labelled by n = 1, 2, . . . , N, with interactions mediated by linear coupling between ball displacements (representing the vibra-tional coupling through a crossbar). In this case we have a network dynamics governed by the equationü where t n i represents the ith impacting event time of the nth node, implicitly determined by u n (t n i ) = ±a. The parameter σ represents a global coupling strength, whilst the specific influence of node m on node n is determined by the value w nm . The network structure is effectively encoded by a matrix with elements w nm . The model equations (4.1-4.2), or variants thereof, also arise naturally when considering mechanical vibro-impact chain systems (Gendelman & Manevitch, 2008;Gendelman, 2013;Perchikov & Gendelman, 2015;Grinberg & Gendelman, 2016). It is convenient to write the network model in first order form by introducing the state vector x n = (u n , v n ) T (where v n =u n ) so thaṫ Here H : R 2 → R 2 describes the form of interaction between the components of nodes and for the (linear) case considered here it is given simply by H(u, v) = (0, u) T . The vector field F : R 2 → R 2 is the single-node dynamics prescribed by F(x n ) = Ax n + f e n , with A as in (2.5) and f e n = (0, f ) T sgn(v n ).
From the form of coupling in (4.3) it is apparent that if x m = x n for all pairs (m, n) then the coupling has no effect and the network reduces to an uncoupled system of individual Franklin bells. Thus, if an individual bell can oscillate then a synchronous network state defined by the N − 1 constraints x 1 (t) = x 2 (t) = · · · = x N (t) = s(t) is guaranteed to exist, where s(t) = (u(t), v(t)) T is the periodic orbit of an isolated node. The techniques for constructing this are precisely those of Section 2.1. The network impact times are also inherited directly from the periodic orbit of an isolated node so that To determine the stability of the synchronous network state it is first convenient to rewrite (4.3) using the graph Laplacian G with components G nm = −w nm + δ nm N k=1 w nk . The network dynamics between impacts then takes the succinct forṁ (4.5) If we now consider a small perturbation to the synchronous orbit by writing x n (t) = s(t) + δx n (t) then we obtain the variational equation where DF(s(t)) and DH(s(t)) are the Jacobian matrices of F(s(t)) and H(s(t)). The PWL nature of the network model means that these can be explicitly calculated as DF(s(t)) = A and DH(s(t)) = 0 0 1 0 . (4.7) If we introduce the vector U = (δx 1 , δx 2 , . . . , δx N ) ∈ R 2N and use the tensor product ⊗, the variational problem (4.6) can be written as Assuming that G is diagonalizable then we can introduce a matrix P formed from the normalized eigenvectors of G such that In each block we have a 2 × 2 linear dynamical system parametrized by the eigenvalues of the graph where ξ l (t) ∈ C 2 . The evolution of the perturbations through the impacting manifolds, where u n = ±a can be obtained using the same approach as in Section 2.1. At the network level this yields given by (2.12). In the transformed coordinates we have that . Thus, we may solve the set of Floquet equations given by (4.10) using the same techniques as deployed for a single node in Section 2.1. This yields ξ l (Δ) = Q(l)ξ l (0), l = 1, . . . , N, where (4.11) One of the eigenvalues of the graph Laplacian is zero (which we fix with the choice λ 1 = 0), with corresponding eigenvector (1, 1, . . . , 1)/ √ N tangential to the periodic orbit. Thus, the synchronous state will be stable if all the other eigenvalues of Q(l), l = 2, . . . , N lie within the unit disc, and the periodic orbit of an isolated node is stable. Since this argument is valid for an arbitrary graph Laplacian it is useful to consider a Floquet problem obtained from (4.11) under the replacement σ λ l → η ∈ C, so that (4.12) If we denote an eigenvalue of Q(η) by q(η) then the MSF is the largest number in the set Re (log q(η))/Δ, and the synchronous state is stable if the MSF is negative at all the points, where η = σ λ l ≡ η l . Thus, the MSF can be computed independently of the network choice and then used to assess the stability of the synchronous state in an arbitrary network, simply by determining where the spectrum of the graph Laplacian lies in relation to the MSF. If any of the eigenvalues of the graph Laplacian lie in the region where the MSF is positive then synchrony is unstable.
We note that it is also natural to consider the stability of the synchronous state in networks of identically coupled limit cycle oscillators using weakly coupled phase oscillator theory. Doing so would give rise to a Kuramoto type network model. The Jacobian determining the stability of the synchronous state would have eigenvalues −σ H (0)λ l , where H (t) is a derived Δ-periodic phase interaction function, and see (Ashwin et al., 2016) for a further discussion. In this case, the stability of the synchronous state is independent of the strength of interaction (though will depend on the network graph Laplacian and the sign of σ ). Thus, it cannot be used to predict any strong coupling instabilities, whereas the MSF approach can.

Master stability function for a smoothed Franklin bell network
The argument above shows how the MSF, originally developed for the study of smooth systems, can be modified for the use of non-smooth systems using saltation operators. We can also sidestep the need to use saltation operators using the smoothing technique described in Section 3. In essence, this leads to the replacement of the saltation operators by propagators exp[A L,R δt] for some fixed small δt with A L,R determining the augmented dynamics in the region, where |u n | > a. In this case the MSF can be constructed in an almost identical fashion to that of Section 4.1 under the replacement of (4.12) by A comparison of the MSF for the non-smooth and smoothed model is shown in Fig. 7. The white region indicates where the MSF is negative. It can be seen that as δt is chosen to be smaller and smaller that the agreement between the two MSFs becomes closer and closer.

Examples
In the following we will illustrate the above concepts with two kinds of network. We begin with an undirected ring network, for which the symmetric coupling strength between nodes n and m is given by where c n , n = 0, 1, . . . , N is non-zero. Due to symmetry, all eigenvalues of the matrix W (with components [W] nm = w nm ) are real. For a network of regular springs, i.e. when c n > 0 for 0 n N, this entails that the eigenvalues of the graph Laplacian are all larger or equal to zero. As a consequence, the synchronous network state is linearly stable since the MSF is negative for all arguments on the positive real half-line, see Fig. 7(D). In Fig. 8(A), we superimpose the η l for a network of 15 nodes, where c n = 1 if n is odd and c n = 0.1 if n is even, while Fig. 8(B) shows results from direct numerical simulations. As expected, the synchronous network state is stable.
It is now instructive to change one of the spring constants c n in the above network to a negative value, which represents a repulsive spring. When we choose c 2 = −0.1, we obtain the results depicted in Fig. 9. In this case, the MSF is negative for one of the η l , say η k , indicating that the synchronous   network state is unstable. Indeed, numerical simulations clearly show a modulation of the values for u n across the network, see Fig. 9(B).
To predict the shape of the emergent network pattern, we can make use of the eigenvector that corresponds to the eigenvalue associated with η k . As Fig. 10 illustrates, the eigenvector resembles very closely the observed values of u n . Note how well the eigenvector captures the large peak and the small oscillations of the network state.
We can now move away from real eigenvalues of W by considering the directed network shown in Fig. 11. The coupling strengths are given by  where μ is a real number, resulting in complex eigenvalues for W. If we choose μ −1 = 2.1 the MSF is negative at the corresponding values of η l , indicating that the synchronous network state is stable (see Fig. 12(A)). Numerical simulations plotted in Fig. 12(B) confirm the results from the linear stability analysis. Here, we plot the time evolution of the v component of all three nodes, i.e. v 1 , v 2 and v 3 . Because of synchrony, the curves overlap and we can see only one trajectory.
When we change μ −1 to 1.9 we obtain a pair of complex conjugates η l that lie in the green region in Fig. 13(A). Here, the MSF is positive, which means that the synchronous network state is unstable. This can also be seen in Fig. 13(B), where we plot trajectories for v 1 , v 2 and v 3 from numerical simulations. In contrast to Fig. 12(B), all three trajectories can be clearly distinguished. Note that the emergent pattern can be predicted from the real eigenvector that is associated with the pair of complex eigenvalues η l . This is an example of a strong coupling instability.

Discussion
Since their inception, Franklin bells have provided the blueprint for numerous electro-mechanical impact oscillators (Asano, 1975;Isacsson et al., 1998;Disna Jayampathi Karunanayake & Hoshino, 2010; Knutson et al., 2007). In its original incarnation, a Franklin bell consisted of two bells between which a metal ball was suspended. The regular chime of a Franklin bell results from the periodic motion of the metal ball between the two bells. Upon impact, the metal ball loses some of its energy, which is captured by a restitution coefficient k < 1. In this regime, only one periodic orbit of the underlying dynamical system (2.3) and (2.4) exists, which is linearly stable (see Fig. 3). As we increase k past one, an unstable solution emerges, which eventually collides with the stable periodic orbit in a saddle node bifurcation. A restitution coefficient larger than one corresponds to an active impacting surface where energy is transferred into the metal ball instead of it being lost from it (Vakakis, 2001;Gendelman, 2006;Pring & Budd, 2011).
For constructing periodic solutions, the non-smooth character of the governing equations does not pose any difficulties. Indeed, we can construct solutions between impacts and then glue them together. Since the system in (2.3) and (2.4) is PWL, solutions are given explicitly in terms of exponential functions. To assess linear stability, we use saltation matrices to propagate perturbations through the impacting manifolds. One could now argue that at a microscopic scale, the dynamics of Franklin bells are actually smooth and the non-smooth character only emerges due to the coarse-grained use of a restitution coefficient. Motivated by this notion, we developed a novel smoothing technique, which is based on supplementing the original dynamical system with two additional parts that describe the dynamics for u > a and u < −a. In each region, we prescribe a linear dynamical system whose coefficients are uniquely determined by demanding that the new pieces of the orbit connect to the existing parts in a C 1 fashion and satisfy the restitution condition. What we need to prescribe, however, is the time-of-flight δt in these two regions. In other words, once we impose a time-of-flight, all coefficients are fixed. The advantage of this approach is that we can explore how the smooth dynamical system approaches the non-smooth one by reducing δt. As Fig. 5 illustrates, letting δt go to zero reduces the propagator in the regions u < −a and u > a to the saltation matrices of the non-smooth system, highlighting the consistency of our new technique.
The above discussion about smooth versus non-smooth representations ties into the discourse on hard impact modelling (particle exposed to a rigid constraint) and soft impact modelling (particle exposed to an elastic constraint) (Blazejczyk-Okolewska et al., 2010). Inelastic models are based on Newton's law of impact and use two main assumptions; (i) the interaction time with the rigid constraint is infinitely small and (ii) that energy dissipation is characterized by a constant restitution coefficient. Explicitly, a restitution coefficient refers to the ratio of post-and pre-impact velocities. In this modelling regime the rigid particle collides with the stiff constraint and none of them are deformed during the collision. On the other hand, soft modelling assumes a finite non-zero contact time and a penetration of the constraint by the colliding body. In this modelling philosophy, the hard impacting constraint is replaced with a spring-damper support or cushioned as it is common in engineering. Elastic impact modelling can be used to analyse different types of spring-damper support systems, which can be either linear or nonlinear (Jiang & Wiercigroch, 2016;Ma et al., 2006Ma et al., , 2008Serweta et al., 2014;Rebouças et al., 2019). Interpolating between these two scenarios is the case in which an elastic body impacts on a rigid surface, which again leads to a nonzero interaction time (Cross, 1999).
While we use the time-of-flight δt to control the transition from smooth to non-smooth dynamics, applications in engineering typically adjust the parameters of the spring-damper system. Naturally, these two approaches are equivalent. Shaw & Holmes (1983) observed that as the stiffness of the cushioned constraint approaches infinity, collision times go to zero and the system becomes an inelastic impact oscillator. Further evidence for this equivalence is provided in Ing et al. (2006Ing et al. ( , 2008. The findings in Blazejczyk-Okolewska et al. (2010) and Jiang et al. (2017) also demonstrate that the smooth system approaches the non-smooth one for large spring-damper stiffnesses. In addition, these studies show how the dynamics of the two systems diverges for softer spring-dampers. This is attributed to the growing influence of external forces, such as gravity, and is consistent with the idea of larger impact times, since only then have these external forces sufficient time to interact.
Having investigated a single Franklin bell, we next turned to networks of N Franklin bells with arbitrary topology. Crucially, each node of the network corresponds to a Franklin bell, and nodes are coupled via springs. Our work contrasts that in Shiroky & Gendelman (2016), where in a linear chain only the central node was a Franklin bell, whilst the remaining nodes were classical non-impacting pendula. Our interest was in the linear stability of the synchronous network state. The existence of synchrony is guaranteed due to the linear coupling between Franklin bells. For linear stability, we employed the MSF approach (Pecora & Carroll, 1998), which reduces the complexity of the linear stability analysis from investigating a 2N-dimensional system of coupled equations to N twodimensional systems. As Fig. 7 illustrates, the MSF for the non-smooth model is well approximated by the one for the smoothed dynamics. However, as we make the time-of-flight δt in the additional regions |u| > a larger, the topology of the MSF changes. A new bubble emerges around the origin, and the extended white region of stability shifts to the right, cf. Figs 7(A) and 7(D). For a ring network of standard springs, i.e. with positive spring constants, the MSF predicts that synchrony is stable, which is confirmed by direct numerical simulations (Fig. 8). By changing the spring constant of one of the springs in the network to a repulsive value, one η l crosses into the green region where the MSF is positive, indicating that the synchronous network state is unstable (Fig. 9). This highlights the fact that subtle changes to the network parameters can have drastic consequences for the network dynamics. Close to the onset of instability, only one η l crosses into the region where the MSF is positive. In this case, the eigenvector associated with the corresponding eigenvalue provides a good estimator for the emergent network state as illustrated by Fig. 10. For the examples above, all eigenvalues of the connectivity matrix W are real. By changing the topology of the network, the eigenvalues of the graph Laplacian may also become complex. Again, the MSF predicts correctly the linear stability of the synchronous state, see Figs 12 and 13.
While we focussed our analysis on Franklin bells, the present study more generally furthers our understanding of networks comprised of nodes with non-smooth dynamics. To date, discontinuous and non-smooth dynamical systems have mostly been studied in isolation. Yet, networks are ubiquitous across engineering and the natural and social sciences. It is therefore desirable to expand our toolbox from individual to interacting non-smooth dynamical systems, as recently advocated by Coraggio et al. (2019) for piecewise-smooth systems, with applications in seismology and load balancing in power grids. As we have illustrated, concepts such as saltation matrices, which are useful at the node level, carry over to the network level and expand the applicability of central techniques for smooth dynamical systems, such as the MSF, to non-smooth systems. A possible extension of this work can be achieved by adding time delays (Steur et al., 2014). Moreover, the techniques used here could be adopted to vibro-impact energy harvesting systems (Yurchenko et al., 2017;Afsharfard, 2018) to test the efficiency at the network level, and our smoothing method could be useful for the investigation of new materials such as elastic support for fenders (Sitnikova et al., 2008(Sitnikova et al., , 2010.

Funding
Engineering and Physical Sciences Research Council (grant number EP/P007031/1).