Breathers and other time-periodic solutions in an array of cantilevers decorated with magnets

1 Department of Mathematics, Bowdoin College, Brunswick, Maine 04011, USA 2 Department of Mechanical and Process Engineering (D-MAVT), Swiss Federal Institute of Technology (ETH), 8092 Zürich, Switzerland 3 Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA, USA 4 Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA 01003-4515, USA

too, the relevant model is a mixed FPUT-KG equation. The microcantilever experiments of [21], however, were conducted in a vacuum to eliminate damping from the model considerations, and thus Hamiltonian dynamics were studied. For our experimental set-up, the Hamiltonian variant is a possible approximate model, however the effects of damping are non-negligible due to the scale of the system. Therefore, the focus of this article will be on a damped-driven, mixed FPUT-KG equation. In particular, time-periodic excitations, including breathers, will be studied. In this respect, our paper can be thought of as a damped-driven extension of the paper [20].
The present paper is organized as follows: In Section 2 we describe the experimental set-up which is modeled in Section 3. We explore the linear problem and outline the numerical methods for the computation of periodic orbits in Section 4. After validating the model in Section 5 via direct comparison with experimental results, we present a number of numerical results in Section 6 concerning the existence, stability and bifurcation structure of time-periodic solutions in various set-ups, including homogeneous ones under various excitations and defect ones. Section 7 concludes the paper.

Experimental set-up
The array of cantilevers we consider is fixed to a rigid aluminum slab [ Figure 1 (center)]. The boundaries of the array are also rigid aluminum slabs with magnets at their tops. Using 3D printed holders, we also embed three magnets within the tip of each cantilevers (right, 6 mm × 4 mm × 2 mm from supermagnete.ch) and orient them along the vertical axis such that the cantilevers repel each other. A phospher bronze metal sheet (M3) from Hassler-Profile is used to manufacture the cantilevers. The base of the cantilever array is guided to move freely along the direction of the first flexural resonance mode of the cantilever.
We use a shaker (Brüel & Kjaer Type-4810) driven by an audio amplifier (Topping TP22) to excite the base of the cantilever array. A lock-in amplifier is used to generate the excitation signal and filter the measured velocity of the base (bottom, Zurich Instruments HF2LI). We acquire the velocity of the base using a laser Doppler vibrometer (Polytec OFV-505, Polytec OFV-5000). The positions of the cantilever tips are acquired at 400 fps using a high speed camera (Vision Research Phantom v1610). Subsequently, we use digital image correlation (DIC) to calculate sub-pixel-resolution positions of the cantilever tips (VIC 2D). Additionally, we use an oscilloscope (Tektronix DPO3014) to acquire the time-domain signals.
To characterize the system we start by considering a single cantilever with no magnets on the system boundaries. We measure the geometric parameters of the cantilever using a caliper and the tip mass using a scale (Table 1). By measuring the dynamic response of the single cantilevers we can experimentally characterize the resonance frequency and the damping characteristics of the cantilevers as shown in Figure 2. The damping coefficient γ is estimated from the damping ratio ξ, which is estimated using where ∂φ ∂ f is the slope of the phase with respect to the swept frequency at the resonance frequency f 0 . Using the tabulated material properties of the phosphor bronze (Table 1) and the other geometric parameters, we predict a resonance frequency of 53.66 Hz (using the methods outlined in Section 3). Experimentally, we found a median resonance frequency of 50.41 Hz. We believe the 6% frequency shift is due to the added flexibility from the fixation of the cantilever that uses small screw connections. To optimize the model, we subsequently estimate the stiffness of the analytical model using the measured resonance frequency and the effective mass. . Probability distribution of the cantilever resonance frequency (left), the damping coefficient (center) and the effective mass of the cantilever model (right). We use the median values for Table 1. We then characterize the force deflection relationship of the magnets using an Instron E3000. Figure 3 shows the measured data and an approximation of the force deflection relationship using where A, p are the fitting parameters. The distance d can be separated into two parts, d 0 (the initial distance) and ∆d (additional deflection). See Table 1. We estimate the inter-magnet distance d 0 by measuring the distance between the outermost magnets and averaging this distance after removing one magnet-width (Table 1). 5 10 15 Figure 3. Characterization of the magnetic interaction forces: The black dots and the red crosses are experimental measurements. We exclude the noisy data (red crosses) for the estimation of the fitting function (blue curve). The experimental data is well described by the fitting function which follows the centerline of the measured force.

Model derivation
We first consider a single beam that is cantilevered to a movable base. If the displacement of the base is y(t) and the displacement of the tip of the beam is x(t) then the lumped-parameter model for the dynamics of the beam tip is given by [22], where the dot stands for differentiation with respect to time t and m, γ, κ are the so-called effective mass, damping coefficient, and stiffness, respectively. This model is derived by considering only the fundamental mode of vibration from the Euler-Bernoulli model of a clamped beam. Through this derivation, one finds the effective stiffness to be where Y I is the beam's bending stiffness and L is the length of the beam. If a tip mass is present (which is the case in our system), the effective mass is found by expressing the total kinetic energy of the beam in terms of the velocity at the tip through the Rayleigh quotient [22], wherem is the mass per unit length of the beam and M t is the mass of the tip. The effective damping coefficient is where ξ is the effective damping ratio. Since we will eventually consider an array of beams that are all attached to the same base, it is useful to introduce the relative displacement z(t) = x(t) − y(t), which represents the position of the tip of the beam relative to the base, where the external driving force is given in terms of the base displacement via the expression −mÿ.
We now consider the effect of adding magnets to the tip of the cantilever. The repelling force of two magnets whose center-to-center separation distance is d 0 is given by F = Ad p 0 where A and p are determined empirically, see Section 2. The experimentally measured values for all quantities and the effective parameter values are given in Table 1. If the beams decorated with magnets are separated by a distance d 0 on the base, and we let the quantity u n (t) represent the displacement of the n − th beam tip relative to its location on the base, then we arrive at the model equation of interest, where n = 1, . . . , N and N is the number of cantilevers. The term βu 3 n is a phenomenological term added to account for the possibility of nonlinearity in the beam [22] for the numerical considerations later in the paper. It is assumed that the base is excited harmonically with the driving frequency f d , e.g., y(t) = y 0 cos(2π f d t) where y 0 is the amplitude of the excitation. Thus, the external force coefficient is given by F ext = −y 0 (2π f d ) 2 m. If a base excitation is introduced, the walls are assumed to be immovable so that u 0 (t) = u N+1 (t) = 0. This is the situation for the experimental set-up. We also will explore numerically the possibility of boundary excitation with a non-moving base, in which case u 0 (t) = α cos(2π f d t), u N+1 (t) = 0 and F ext = 0.

Linear analysis and computation of periodic orbits and their stability
The Hamiltonian lattice corresponds to β = y 0 = γ = 0. To consider linearized dynamics (with p 1) we require small amplitudes, i.e., |u n −u n+1 | d 0 1, so that the magnetic coupling force can be linearized. where the dispersion relationship is given by where K 2 = pAd p−1 0 < 0. In the case of N → ∞, i.e., an infinite lattice, k ∈ [0, π], and for a lattice consisting of N (finite) nodes with zero boundary conditions k = k n = nπ N+1 with n = 1, 2, . . . , N. The dispersion curve in terms of frequency f = ω/2π is shown in Figure 4 (left). The lower and upper cutoff frequencies, f = ω(0)/2π ≈ 50.6 Hz and f u = ω(π)/2π ≈ 126.8 Hz respectively, are shown as dashed lines in Figure 4.
In the damped-driven linear system, namely when β = 0, |u n −u n+1 | d 0 1 and y 0 0 γ timeperiodic steady-states of the form u n (t) = a n cos(2π f d t + φ n ) can be found explicitly, where a n and φ n are the amplitude and phase of the nth node, respectively. An example plot of a 5 for f d ∈ [0, 140] and F ext = 1.54 · 10 −6 N is shown in Figure 4 (right). The critical points correspond to the N = 21 linear resonant frequencies. An example linear resonance at f ≈ 75.2 Hz is shown as a red dot in Figure 4 (right).
For the remainder of the article, we will be interested in time-periodic orbits of the nonlinear Eq. (3.3) (i.e., where the conditions β = 0 and |u n −u n+1 | d 0 1 do not necessarily hold). In this case, exact analytical solutions are not available (although they can be approximated via long wavelength approximations [23]). Thus, we will turn to numerical computations. We compute time-periodic orbits of Eq.
Note that the solution frequency and drive frequency are both f d by construction. An initial guess for the Newton iterations can be provided by the linear steady-state solutions u n (t) = a n cos(2π f d t + φ n ). To investigate the dynamical stability of the obtained states, a Floquet analysis is used to compute the multipliers associated with the solutions. The Floquet multipliers for a solution are obtained by computing the eigenvalues of the monodromy matrix (which is V(T b ) upon convergence of the Newton scheme). If a solution has all Floquet multipliers on the unit circle, the solution is called (spectrally) stable. An instability that is result of a multiplier on the (positive) real line is called a real (exponential in nature) instability. However, there can also be oscillatory instabilities, which correspond to complex-conjugate pairs of Floquet multipliers lying outside the unit circle (in the complex plane). In the bifurcation diagrams that are to follow in this paper, solid blue segments will correspond to stable parametric regions while the dashed red segments correspond to real unstable regions and green dashed-dotted segments correspond to oscillatorily unstable regions. The bifurcation diagrams are obtained using a pseudo-arclength continuation algorithm.  We now present experimental data that demonstrate that Eq. (3.3) is a reasonable model for the cantilever array decorated with magnets.

Model validation
In the experiment, we perform amplitude sweeps at a fixed frequency (in particular at f d = 75.2 Hz, which is close to a resonant frequency, see the red dot in Figure 4). In particular, we let the system settle to steady-state for a fixed frequency (and the steady-state amplitude is recorded) and then the drive amplitude is changed slightly and we let the system once again reach steady-state, etc. The experimental results are shown in Figure 5 (right), where the upward sweep is shown as square markers and the downward sweep is shown as triangles. The hysteretic behavior identified in the experiments reveals that the system is bistable for various values of the base excitation amplitude (see e.g., the dashed line at 7µm of Figure 5 (right)). To validate our model, we performed a similar procedure theoretically. In particular, we carried out a pseudo-arclength continuation of time-periodic solutions of Eq. (3.3) with respect to the amplitude y 0 and the drive frequency fixed to f d = 75.2 Hz (other parameter values are given in Table 1). However, no region of bistability is found, suggesting a slight mistuning of our system parameter values. To investigate this further, a pseudo-arclength continuation with respect to the drive frequency f d was performed for Eq. (3.3) with the base excitation amplitude fixed to y 0 = 7µm, see Figure 5 (left). Note the pseudo-arclength continuation allows us to follow theoretically not only the stable branches (solid blue lines) of the bistable response, but also the unstable ones in between (red dashed lines). It is clear from Figure 5 (left) that the region of bistability is for drive frequencies larger than f d = 75.2 Hz. For this reason, we revisited the drive amplitude continuation, but for frequencies slightly larger than 75.2 Hz. In particular, we conducted a pseudoarclength continuation of time-periodic solutions of Eq.

Numerical exploration of the system
Since we demonstrated in the previous section that Eq. (3.3) is a good model for the array of cantilevers decorated with magnets, we will proceed to carry out a more detailed numerical exploration of the system. In the subsequent section, comparisons of boundary and base excitation will be made. In Section 6.2 breather solutions of a lattice with a defect are explored.

Comparing boundary and base excitation
To draw comparisons to systems that are typically driven on the boundary, such as the granular chain, we will now conduct numerical continuations with F ext = u N+1 = 0 and u 0 = α cos(2π f d t). The parameter values for the subsequent comparison of boundary and base excitation are shown in Table 2. We chose different parameter values than in Section 5 to demonstrate how the solutions of the base excitation problem can vary with a different parameter set.
A linear analysis of the resonant peaks will reveal that the most dominant mode in the base excitation problem corresponds to in-phase motion (associated with the lowest frequency response peak), as expected. The modal contributions in the boundary excitation problem are more distributed. This is also reflected in the resonant peaks of the nonlinear system which were found by performing a pseudo-arclength continuation in the drive frequency f d , see Figure 6. For the base excitation, the amplitude is fixed to F ext = 0.028 N (recall that F ext = −y 0 (2π f d ) 2 m). To make a comparison of the same magnitude in the boundary excitation problem, the amplitude is fixed to α = 46.368 µm, which corresponds to an approximate maximum force contribution of pA(d 0 ) p−1 α ≈ 0.022 N. In addition to how the magnitude of the resonant peaks differs between the two excitation problems, the boundary    excitation problem differs in its stability properties. In the base excitation problem, no regions of oscillatory instability were identified (as seen in Figure 6 (top), but also in other similar computations we have made). This is in contrast to the boundary excitation problem, where oscillatory instabilities appear generically (as seen in Figure 6 (bottom), but also in other similar computations we have made). Another important difference is that the instability regions invade (some of) the higher amplitude branches of the resonant peaks in the boundary excitation problem (in Figure 6 (bottom) oscillatory instabilities are present for higher frequency peaks). Thus, while the distribution of the modes in the boundary excitation problem is potentially more appealing from a potential application perspective (for example in the energy harvesting context [22]), this problem is more unstable than the base excitation problem. This latter aspect is even more apparent when inspecting the bifurcation diagram for the amplitude continuation (with f d fixed), see Figure 7. In the base excitation problem the higher amplitude branch is stable, whereas in the boundary excitation problem the higher amplitude branch becomes unstable (via an oscillatory instability) at some critical amplitude (in the case of the right panel of Figure 7 at α ≈ 60) and does not regain stability. Example solution profiles (corresponding to the letter labels of Figure 7) and their Floquet multiplier spectra for the base and boundary excitation problems are given in Figures 8 and 9, respectively.  Figure 6. The letter labels in both panels correspond to the solution profiles shown in Figures 8 and 9. In the simulations shown here, there are N = 8 nodes.

Breathers in a defect lattice
We will now turn our attention to breather solutions of this system, namely ones that are timeperiodic and localized in space. To introduce a localized mode into the system, a light mass defect is added to the center of the chain. In particular, the effective mass of all the particles is the same as before (namely the effective mass is 3.59 g) however the center particle has an effective mass of 1.795 g.
Since the breather solutions we are seeking will be a result of this mass defect, the breathers can also be thought of as nonlinear localized modes [24]. By the term "breathers", we will refer to exponentially localized in space, periodic in time solutions more generally, without necessitating that the underlying lattice system be homogeneous in space. While breathers that are the result of modulation instability (MI) of plane waves from the bottom of the spectral band are possible in the Hamiltonian system (also called intrinsic localized modes), our system is highly damped, and thus obtaining breathers through MI was not possible. Other than this mass defect, and the lattice size which we take to be N = 21 in this section, all parameters are the same as those in Table 2 and we only consider base excitation. The nonlinear resonant peaks of this system are shown in Figure 10. Similar to the top panel of Figure 6, the peak corresponding to the in-phase mode is fairly dominant. However, there is now a resonant peak in the gap, which corresponds to the localized defect mode. A zoom of the defect mode resonant peak is shown in the left panel of Figure 11. With the drive frequency fixed to f d = 174 Hz (see the dashed line in the left panel of Figure 11) an amplitude continuation is also performed, see the right panel of Figure 11. Note, in [25] a similar situation was explored in the context of a granular chain with a single light mass defect that is driven on one boundary. There the higher amplitude solutions in the drive amplitude continuation are unstable (similar to the right panel of Figure 7). At the drive amplitudes where the breathers are oscillatorily unstable, quasi-periodic solutions are approached. In our base excited system, the high amplitude breathers are stable, and thus, the system is not driven to quasi-periodicity. Example breather profiles (corresponding to the letter labels of the right panel of Figure 11) and their Floquet multiplier spectra are given in Figure 12. Note the asymmetry of the spatial profiles of the solutions shown in Figure 12. This is a consequence of the asymmetric nature of the interaction potential. This can be seen through a Taylor expansion of the magnetic force, leading to both quadratic and cubic force terms. In nonlinear lattices with cubic terms only, the potential is symmetric, and hence the profiles are also symmetric [26].
The final numerical experiment considers non-zero values of the onsite nonlinear force coefficient β. Bistability of the defect breather is possible for zero, negative, and positive values of β, see Figure 13. Example breather profiles (corresponding to the letter labels of Figure 13) and their Floquet multiplier spectra are given in Figure 14. Note the bifurcation values (i.e., the turning points in Figure 13) occur for very large values of β. This is due to the fact that the beam deflections are somewhat small, and hence a large value of β is required for the nonlinear term to have an effect on the system. A more physically realistic scenario where the beam nonlinearity would be observed is in the case of large deflections. This would require the equilibrium distance d 0 to be larger (in order to allow "room" for    Here we see that if we increase the dimensional d 0 by, e.g., a factor 15, then according to the above, the coefficient for the onsite cubic nonlinearity is 8 orders of magnitude larger (note the amplitude of the external force increases about 4 orders of magnitude). Thus, while the bifurcation values in Figure 13 seem very large, they have a reasonable physical interpretation, if e.g., considering large deflections. Although such large deflections may be hard to realize with the cantilevers considered in our experiment, other similar systems, such as the one the studied in [27], would allow such large deflections.

Conclusion
We have introduced a system consisting of an array of cantilevers that are coupled via magnetic links and demonstrated that the damped-driven mixed FPTU-KG model provides a reasonably good description of the dynamics observed in the experiments. This article focused on the topic of timeperiodic solutions (and their breather variants), both for the case of base and boundary excitation, and also in the presence of a light mass defect. The bistable response of the system and the possibility of hysteresis were featured both in our numerical and in our experimental results, in qualitative agreement between the two.
Nevertheless, there are many other avenues for future research that this system inspires. One includes revisiting the Hamiltonian variant of the problem to study how the stability of Hamiltonian breathers is affected as the strength of the nonlinear parameter β varies and contrast this to the case where the intersite nonlinearity is a polynomial (as in [20]) as opposed to a power-law, as in our paper. Higher dimensional analogs of this system, in e.g., hexagonal or square configurations, would also be worth exploring. Finally, the robustness of this system for the purposes of applications, such as energy harvesting, would be another important research direction.