Nonlinear Localized Modes in Two-Dimensional Hexagonally-Packed Magnetic Lattices
C. Chong, Yifan Wang, Donovan Marechal, E. G. Charalampidis, Miguel Moleron, Alejandro J. Martinez, Mason A. Porter, P. G. Kevrekidis, Chiara Daraio
NNonlinear Localized Modes in Two-Dimensional Hexagonally-Packed Magnetic Lattices
C. Chong*, Yifan Wang, Donovan Mar´echal, E. G. Charalampidis, Miguel Moler´on, Alejandro J. Mart´ınez,
5, 6
Mason A. Porter, P. G. Kevrekidis, and Chiara Daraio Department of Mathematics, Bowdoin College, Brunswick, ME 04011, USA Division of Engineering and Applied Science California Institute of Technology Pasadena, CA 91125, USA Mathematics Department, California Polytechnic State University, San Luis Obispo, CA 93407-0403, USA Institute of Geophysics, Department of Earth Sciences, ETH Zurich, 8092 Zurich, Switzerland Computational Biology Laboratory, Fundaci´on Ciencia & Vida, Santiago, 7780272, Chile Universidad San Sebastian, Santiago, 7510156, Chile Department of Mathematics, University of California, Los Angeles, CA 90095, USA Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA, 01003, USA
We conduct an extensive study of nonlinear localized modes (NLMs), which are temporally periodic andspatially localized structures, in a two-dimensional array of repelling magnets. In our experiments, we arrange alattice in a hexagonal configuration with a light-mass defect, and we harmonically drive the center of the chainwith a tunable excitation frequency, amplitude, and angle. We use a damped, driven variant of a vector Fermi–Pasta–Ulam–Tsingou lattice to model our experimental setup. Despite the idealized nature of the model, weobtain good qualitative agreement between theory and experiments for a variety of dynamical behaviors. We findthat the spatial decay is direction-dependent and that drive amplitudes along fundamental displacement axes leadto nonlinear resonant peaks in frequency continuations that are similar to those that occur in one-dimensionaldamped, driven lattices. However, driving along other directions leads to the creation of asymmetric NLMsthat bifurcate from the main solution branch, which consists of symmetric NLMs. When we vary the driveamplitude, we observe such behavior both in our experiments and in our simulations. We also demonstrate thatsolutions that appear to be time-quasi-periodic bifurcate from the branch of symmetric time-periodic NLMs.
I. INTRODUCTION
Discrete breathers are spatially localized, time-periodic solutions of nonlinear lattice differential equations. They have beenstudied in a host of scientific areas, including optical waveguide arrays and photorefractive crystals [1], Josephson-junctionladders [2, 3], layered antiferromagnetic crystals [4, 5], halide-bridged transition-metal complexes [6], dynamical models of theDNA double strand [7], molecular lattices [8], Bose–Einstein condensates in optical lattices [9], and many others.Most of the immense volume of work — now spanning more than three decades — on discrete breathers has been in one-dimensional (1D) lattices [10–12]. Most relevant to the present article is research on discrete breathers in Fermi–Pasta–Ulam–Tsingou (FPUT) lattices, which have nonlinear inter-site coupling [13, 14]. Moreover, FPUT-like lattices with power-law poten-tials have been used to model a variety of mechanical systems, such as granular crystals [15–18] and (more recently) magneticlattices [19–21].There have also been some studies of breathers in two-dimensional (2D) lattices, although there are many fewer such studiesthan of 1D lattices. Example physical settings in 2D include crystal lattices [12, 22], electric circuits [23], and dusty plasmas[24, 25]. Breathers in 2D lattices have been analyzed with both asymptotic methods [26] and numerical methods in bothhomogeneous [27, 28] and heterogeneous media [29, 30]. See [8, 31] for overviews of results about 2D breathers.The 2D setting of the present work is a mechanical system in which each magnet has two displacement fields, which distin-guishes it from many studies of scalar 2D lattices, such as those that describe electrical circuits [23]. Specifically, we examine alattice of repelling magnets that are arranged in a hexagonal configuration. The choice of a hexagonal arrangement is motivatedby the experimental setup, as hexagonal configurations are more robust structurally than other arrangements (such as squareones). At the center of the lattice is a light-mass defect, which introduces a localized defect mode into the spectrum of thelinearization of the system. To excite the system experimentally, we drive the center of the lattice by a force that results from thecurrent that flows along a wire that we suspend above the lattice. We model damping using a dashpot term. Putting everythingtogether, the proposed model for the experimental setup is a damped, driven variant of a vector FPUT lattice.Although a breather is defined as a spatially localized and time-periodic structure, it is useful to label different types ofbreathers. Linear systems with an impurity or a defect (e.g., with a particle of lighter mass than the other particles) have isolatedpoints in their spectra that lie above the spectral edge. These modes are called “defect localized modes” [32]. In the presenceof nonlinearity, breathers can bifurcate from these modes and can exist for frequencies other than the linear defect frequency.Breathers that manifest in this way are called “nonlinear localized modes” (NLMs) [33] and are the subject of the present work.By contrast, breathers that do not manifest via a defect or an impurity are called “intrinsic localized modes” (ILMs). One wayfor ILMs, which we do not investigate in the present paper, to manifest themselves is via a modulation instability of plane waves[10]. In addition to breathers, other kinds of orbits — such as quasiperiodic and chaotic ones — can also occur in nonlinearlattices. For example, such orbits have been identified in strongly nonlinear damped, driven granular chains [34, 35], suggesting a r X i v : . [ n li n . PS ] S e p (a) (b) (c) Displacement (mm) F o r c e ( N ) FIG. 1: (a)
Picture (from a video frame that was used for particle motion tracking) of the experimental setup. The yellow arrow indicates thedirection of the external excitation. (b)
Sketches of the (top) normal and (bottom) defect particles. (c)
Magnetic dipole–dipole interactions inexperiments (open grey circles) and fitted model using Eq. (2) (solid black curve). that such solutions may also be present in damped, driven magnetic lattices. In the present work, we examine such NLM states,their stability, and the modes that arise as a result of instabilities.Our paper proceeds as follows. We present our experimental setup in Sec. II, and we detail the corresponding model equations,linear theory, and numerical methods in Sec. III. We give the main numerical and experimental results in Sec. IV, where weexplore NLM profiles, spatial decay, parameter continuations, and nearly time-quasiperiodic orbits. We conclude and discussfuture challenges in Sec. V.
II. EXPERIMENTAL SETUP
We place a 2D lattice on an air-bearing table to make the magnetic particles (which constitute the nodes of the lattice) levitate.The lattice consists of 127 magnetic particles that are hexagonally packed. We glue 36 of these particles to the boundaries,and 91 of them are free to move [see Fig. 1(a)]. Each particle is a 3D-printed disc with a hole in the center, where we attacha neodymium magnetic cylinder. We glue a thin piece of cover glass at the bottom to make the surface smoother and therebyimprove the levitation of the particles. We build the defect particle, which is located in the center of the lattice, by directlyattaching the magnet on the glass without the 3D-printed structure. This particle has a lighter mass and serves as a defect [seeFig. 1(b)]. The mean mass of a normal disc particle is 138.2 mg ± particles). The defect particle has a mass of 81.6 mg, which corresponds to 58.68% of the normal particle mass.We excite the defect particle using an external magnetic field that we generate using a conductive wire that we place over theparticle at a height of mm. We generate the AC current that flows through the wire from a Lock-In amplifier (SR860 500kHzDSP Lock-in Amplifier), and we amplify it with an audio amplifier (Topping TP22, class D). The equation that describes theforce that the wire exerts on a magnet at distance r from it can be calculated as: F wire ( r ) = Iµ M π h − r ( h + r ) , (1)where h is the height of the wire from the plane of floating discs, I is the wire current, µ = 4 π · − N A − is the magneticpermeability, and M = 7 . · − Am is the magnetic moment of the floating disc. See the appendix for the derivation ofEq. (1). We use harmonic excitations in our experiment, so the current through the wire is I ( t ) = aI sin(2 πf t ) , where f is thedrive frequency (in Hz), a is the drive-voltage amplitude (in Volts), and I = 0 . A V − is the current per unit voltage that wemeasure in the wire.The magnets repel each other. In the ideal situation of a perfect dipole–dipole interaction, the magnetic force between tworepelling magnets is F magnet ( r ) = Ar p , (2)where r is the distance (in meters) between the two center points of the magnets, p = − , and A = 3 µ M / π . AlthoughEq. (2) is reasonable for large separation distances, we obtain better agreement by empirically determining A and p . Becausethe force between two magnetic dipoles is too small to measure directly, we create a pair of plastic plates, with 25 magnetsattached to each plate. We position the plates to align each pair of cylindrical magnets from opposite plates through their radialdirections. We measure the repulsive force as a function of the displacement between these two plates in a materials tester(Instron ElectroPuls E3000). The distance between the magnets on each plate is large enough (specifically, it is 2.5 cm) so thatwe can neglect interactions between magnets that are not aligned. The distance between a magnet on the first plate and thenon-aligned magnets of the other plate is larger than 25 mm. As one can see in Fig. 1(c) the interaction force already tendsto for distances smaller that are significantly smaller than 25 mm. Consequently, the measured force is approximately equalto the sum of the repulsive force of the 25 isolated magnet pairs. We fit the data using Eq. (2), which yields p ≈ − . and A ≈ . · − N/m p [see Fig. 1(c)].We monitor the motion of the central particle using a laser vibrometer (Polytec CLV-2534), and we record the remainder ofthe lattice using a digital camera (Point Grey GS3-U3-41C6C-C) with a frame rate of 90 fps. We analyze the images usingdigital-image-correlation (DIC) software (VIC-2D) to determine each particle’s velocity. We inspect half of the lattice, as thecables that are connected to the driving wire block most of the system’s other half [see Fig. 1(a)]. Due to imperfections at thebottom of the glass disks (e.g., dust, scratches, and so on) and the fact that mass is not distributed evenly on a disk, a few particlesstart to rotate when they are levitated by the air that flows out of the air-bearing table. The image correlation software then losestrack of them. We ignore these rotating discs in our subsequent analysis. To estimate the value of the damping coefficient, weexcite the center particle and record the resulting temporal amplitude decay once we switch off the excitation. We fit the decayto an exponential function, which we then use to calculate the damping coefficient γ ≈ . · − N s/m. The lattice particlesare always in motion with at least small speeds, even in the absence of excitation. This is due to interactions with the air flowfrom the table and to imperfections (e.g., nonaxisymmetric mass distributions) of the particles. We use this motion to estimatethe noise in the system. To evaluate the amount of noise, we record the lattice motion without excitation as a comparison. Wesummarize the values of the parameters in Table I.
TABLE I: Summary of the parameter values in our experimental setup.Description Symbol Value (measured) Description Symbol Value (fitted)Mass of bulk magnet M b . mg Magnetic coefficient A . · − N/m p Defect mass M δ . mg Nonlinearity p − . Equilibrium distance δ . mm Damping coefficient γ . · − N s/mWire height h mm Magnetic moment M . · − A m III. THEORETICAL SETUPA. Model Equations
Our goal is to study NLMs in a 2D hexagonal lattice. In selecting equations to model the system that we described in Sec. II,we seek the simplest possible model that incorporates the ingredients (nonlinearity, discreteness, and dimensionality) that areessential for NLMs and also yield reasonable agreement with experimental data. It is in this spirit that we develop our modelequations. After doing so (at the end of this section), we briefly discuss model simplifications.We consider a hexagonally packed lattice of magnets. We use the lattice basis vectors e = (1 , and e = (1 / , √ / .Let q m,n ( t ) = ( x m,n ( t ) , y m,n ( t )) ∈ R denote the displacement from the static equilibrium of the magnet at position p = δ ( me + ne ) in the plane [see Fig. 2(a)], where δ is the center-to-center distance between two particles at equilibrium. Thelattice indices m and n take the values m, n ∈ {− w, − ( w − , . . . , , . . . , w − , w } , where w is the number of magnets alongan edge of the hexagon. The lattice boundary is given by the hexagon with magnets at positions ( wδ cos( jπ/ , wδ sin( jπ/ ,where j = 0 , , . . . , w − [see Fig. 2(b)]. For our fixed boundary conditions along the edge of the hexagonal boundary, q m,n ( t ) = 0 if | m + n | > w .One can express the distance between the magnet with index ( m, n ) and one of its nearest neighbors in terms of the displace-ments x m,n and y m,n of the magnets from their respective equilibrium positions. Once we determine the distance, we computethe resulting force using Eq. (2). Summing the forces from each of the six nearest neighbors and applying Newton’s second lawleads to the following equations of motion: M m,n ¨ q m,n = − F ( q m +1 ,n − q m,n ) − F ( q m,n +1 − q m,n ) + F − ( q m,n − q m − ,n +1 )+ F ( q m,n − q m − ,n ) + F ( q m,n − q m,n − ) − F − ( q m +1 ,n − − q m,n ) − γ ˙ q m,n + F ext m,n ( t ) . (3) ( m, n ( m + 1 , n m , n ) ( m, n ) ( m + 1 , n )( m , n + 1) ( m, n + 1) mn (a) (b) FIG. 2: (a)
Orientation for our convention of indexing particles in a hexagonal lattice. The m axis and n axis meet at an angle of θ = π/ . (b) A hexagonal lattice with w = 6 magnets along each edge of the lattice. The empty circles and solid center circle represent the locations ofthe magnets in equilibrium. When in equilibrium, the center-to-center distance between any two adjacent magnets is δ . The outer hexagonalboundary is the solid gray hexagon that encloses the lattice. On the boundary, the solid points represent fixed (i.e., immovable) magnets. Thereare w + 1 such magnets along each edge of the boundary. The solid circle represents the defect particle, which has index ( m, n ) = (0 , . The vector functions F j ( q ) = F j ( x, y ) ∈ R have a magnitude of | F j ( x, y ) | = A (cid:20)(cid:113) ( δ cos( θ j ) + x ) + ( δ sin( θ j ) + y ) (cid:21) p , θ j = jπ , j ∈ {− , , } . The mass of the magnet with index ( m, n ) is M m,n . The dashpot γ ˙ q m,n is a phenomenological term that we add to accountfor damping. Using such a term has yielded reasonable agreement with experiments in other, similar lattices [19, 36, 37]. Thequantity F ext m,n is the external force that we apply to the magnet at ( m, n ) . In the present article, we consider excitations via awire that is directly above the center of the lattice. The magnitude of the excitation is given by Eq. (1). Therefore, F ext0 , ( t ) = a sin(2 πf t ) I µ M π cos( φ ) h − x , ( h + x , ) sin( φ ) h − y , ( h + y , ) , (4)where φ is the angle of the excitation and F ext m,n = 0 when m (cid:54) = 0 and n (cid:54) = 0 . In our experiments and in most of our numericalcomputations, the excitation angle is φ = π/ , so we excite only the y -component of the center magnet. We will also exploresome other excitation angles. As we discuss in the appendix, the lattice forces dominate the dynamics. The wire has only a smalleffect on magnets that are beyond the center of the lattice. For example, at equilibrium, the force that is exerted on the centermagnet by the wire is two orders of magnitude larger than the force that the wire exerts on the center magnet’s nearest neighbors.Compare Eq. (1) with r = 0 and r = δ .In our model, we ignore effects beyond nearest-neighbor coupling of the magnetic interactions. It is known that such long-range effects can alter the structure of localized modes. For example, it was shown in [38] that the spatial decay of breathers cantransition from exponential spatial decay to algebraic decay in lattices with algebraically decaying interaction forces (as is thecase in our model) for lattices with sufficiently many sites. More recently, Moler´on et al. [37] studied NLMs in a 1D magneticlattice using a model with long-range interactions. Although the differences between long-range and nearest-neighbor latticesthat were considered in [37] are detectable, they are still small. For example, at equilibrium, the force that is exerted on thecenter magnet by its nearest neighbors is one order of magnitude larger than that exerted by its next-nearest neighbor. (CompareEq. (2) with r = δ and r = 2 δ .) Therefore, to keep our model as simple as possible, we ignore such small long-range effects.In our analysis of experimental data we ignore magnets that are rotating, so our model does not account for rotation. Thisleaves air resistance as the primary source of damping. Given the size of the magnets and velocities that we consider, we employa linear dashpot [39]. We also assume that the magnets stay in a plane. We validate the many assumptions that we have made informulating the model in Eq. (3) via a direct comparison with experimental results in Sec. IV.For the remainder of the manuscript, we fix all parameters of the model (and we summarize them in Table I), except for theexcitation amplitude a , frequency f , and angle φ . We will specify these in our various examples. In all cases, we examinea lattice with a single defect particle in the center and a hexagonal boundary with a length of w = 6 magnets (see Fig. 2).Importantly, we do not fit the parameter values to the reported experimental results. Instead, we determine them beforehandusing the procedures that we detailed in Sec. II. -5 0 5 ℓ -505 k KM Γ (a) -5 0 5 ℓ -505 k (b) M K Γ M F r equen cy [ H z ] (c) FIG. 3: (Color Online) (a)
Contour plot of the bottom dispersion surface. We show the irreducible Brillouin zone as the triangle withmagnets at the points that we mark by M , K , and Γ . (b) Contour plot of top dispersion surface. (c)
Band structure along the edge of theirreducible Brillouin zone [also see the triangle in panel (a)] for the bottom (dashed gray curve) and top (solid black curve) dispersion surfaces.The horizontal dashed curve corresponds to the defect-mode frequency f ≈ . Hz in a finite-dimensional system with a mass ratio of M δ /M b ≈ . , where M δ denotes the mass of the defect magnet and M b denotes the mass of the other (“bulk”) magnets. B. Linear Analysis
We start with the basic linear theory of localized modes for our hexagonal magnetic lattice. We are particularly interestedin modes with frequencies that lie above the cutoff frequency of the pass band. We first derive an analytical expression for thecutoff frequency, which is straightforward for an infinite-dimensional Hamiltonian system (i.e., with all integers m and n , alongwith a = 0 and γ = 0 ). We then numerically estimate the frequency of a linear mode that is associated with the defect forin the finite-dimensional Hamiltonian system. Finally, we compute linear localized modes in the associated finite-dimensionaldamped, driven system.Assuming small strains, such that | q m ± ,n − q m,n | δ (cid:28) , | q m,n ± − q m,n | δ (cid:28) , | q m ± ,n ∓ − q m,n | δ (cid:28) , (5)one can Taylor expand, F j ( q ) ≈ F j ( q ) + D F j ( q ) q , where D F j is the Jacobian matrix of F j . Using this notation, we write the linearized equations of motion: M m,n ¨ q m,n = − D F ( q m +1 ,n + q m − ,n ) − D F ( q m,n +1 + q m,n − ) − D F − ( q m − ,n +1 + q m +1 ,n − ) + 2( D F + D F + D F − ) q m,n , (6)where the entries of the Jacobian matrices are denoted by, D F j = a j b j c j d j , j ∈ {− , , } , with a − = p ˆ δ , b − = 0 , c − = 0 , d − = ˆ δ ,a = 3 + p δ , b = √ p − δ , c = b , d = 1 + 3 p δ ,a = a , b = − b , c = − c , d = d , where ˆ δ ≡ Aδ p − . For a monoatomic system (in which all magnets are identical, such that M m,n = M b ), the linear system hasplane-wave solutions q m,n = q exp (cid:16) i ( km + n k + √ l )) (cid:17) e iωt , q ∈ C , k, (cid:96), ω ∈ R , where the wavenumbers k, (cid:96) and angular frequency ω = ω ( k, (cid:96) ) satisfy the dispersion relationship [ ω ( k, (cid:96) )] = ω a + ω d ± (cid:113) ( ω a + ω d ) − ω a ω d − ω b ω c )2 , (7)where ω α ( k, (cid:96) ) = ( − α − cos( k ) − α cos( k/ √ / l ) − α cos( k/ − √ / l ) + 2( α − + α + α )) /M b and α ∈ { a, b, c, d } . In Fig. 3(a,b), we show contour plots of the two dispersion surfaces from Eq. (7). The dispersion curvesalong the edge of the irreducible Brillouin zone are shown in Fig. 3(c). The cutoff value of the pass band has the wavenumberpair ( k, (cid:96) ) = (0 , π/ , which is where the dispersion curve attains its maximum value. For the parameter values in Table I, thecutoff frequency is f c = ω (0 , π/ + / (2 π ) ≈ . Hz, where ω + corresponds to the top dispersion surface.The presence of the lighter defect introduces a linear mode into the system that is localized in space and oscillates with afrequency above the cutoff frequency of the linear monoatomic system. With the light-mass defect at the center of the lattice,we write M m,n = M δ , n = 0 and m = 0 M b , otherwise , (8)where (0 , is the index of the mass defect with mass M δ , the quantity M b is the mass of a magnet in the “bulk” (i.e., thenon-defect mass), and M δ < M b . We numerically compute the linear modes of the system with a mass defect, and we arethereby able to consider finite lattices. We use a hexagonal boundary with edge length w = 6 magnets [see Fig. 2(b)]. One canembed this lattice into a square matrix of size N × N , where N = 2 w − is the number magnets along the n = 0 line of thelattice. Let X ( t ) be the N × N matrix whose ( m, n ) th entry is x m,n ( t ) , and let Y ( t ) be the N × N matrix whose ( m, n ) thentry is y m,n ( t ) . We enforce the fixed hexagonal boundaries by setting the displacements of magnets with indices ( m, n ) suchthat | m + n | > w to . We define the N × N matrix operators L α through L α Y = α DY + α Y D + α ( E T Y E T + EY E − Y ) , (9)where α ∈ { a, b, c, d } ; the N × N tridiagonal matrix D has entries on the super-diagonals and sub-diagonals, entriesalong the diagonal, and entries everywhere else; E is an N × N matrix with entries along the super-diagonal and entrieseverywhere else; and E T is the transpose of E . With these definitions, Eq. (6) becomes M ◦ ¨ X ( t ) = L a X ( t ) + L b Y ( t ) ,M ◦ ¨ Y ( t ) = L c X ( t ) + L d Y ( t ) , (10)where M is an N × N matrix in which all entries except the (0 , entry (which is equal to M δ ) are equal to M b . The operation ◦ denotes pointwise multiplication (i.e., the Hadamard product). The system (10) has solutions of the form X ( t ) = ˜ Xe iωt and Y ( t ) = ˜ Y e iωt , where ˜ X and ˜ Y are N × N time-independent matrices and − ω M ◦ ˜ XM ◦ ˜ Y = L a L b L c L d ˜ X ˜ Y . (11)One can cast Eq. (11) as a standard eigenvalue problem by letting λ = − ω and unwrapping the ˜ X and ˜ Y matrices intoequivalent row vectors and reshaping the block matrix (with entries given by L α ) into a corresponding N × N matrix. Onecan then numerically solve the resulting eigenvalue problem to obtain N eigenvalues and their corresponding modes. Usingthe values in Table I and w = 6 (which yields N = 11 ), we see that two eigenvalues (each with a frequency of √− λ/ (2 π ) = ω/ (2 π ) = f d ≈ . ) lie above the cutoff frequency f c ≈ . Hz. We plot the value of f d as the horizontal dashed line inFig. 3(c), and we show the two corresponding modes in Fig. 4(a,b). Both of these modes are spatially localized, as desired.Now suppose that there is driving and damping. Near the background state, equations (10) yield the following approximatesystem: M ◦ ¨ X ( t ) = L a X ( t ) + L b Y ( t ) − γ ˙ X ( t ) + a A cos( φ ) sin(2 πf t ) , (12) M ◦ ¨ Y ( t ) = L c X ( t ) + L d Y ( t ) − γ ˙ Y ( t ) + a A sin( φ ) sin(2 πf t ) , (13) (a) (b) Frequency [Hz] R M S [ mm / s ] (c) Frequency [Hz] R M S [ mm / s ] (d) FIG. 4: (a,b)
Shape of the two modes with defect frequency f d ≈ . Hz of the linear Hamiltonian system (6) with w = 6 magnets alongeach edge of the boundary. The color intensity at each point ( m, n ) corresponds to (cid:112) ˙ x m,n + ˙ y m,n . (c) RMS of ˙ y , of the linear damped,driven solution (14) as a function of the drive frequency f with excitation amplitude a = 0 . mV and φ = π/ . (d) RMS of the velocity ofthe center particle of the experimental frequency sweep with a = 4 mV and φ = π/ . where A is an N × N matrix that has all entries except for the single nonzero entry A , = I µ M πh . We obtain A by expanding the external drive function F ext near the vanishing displacements and maintaining the leading,non-vanishing term. We can then find solutions of the system (12,13) in the form X ( t ) = ˜ X cos(2 πf t ) + ˜ X sin(2 πf t ) , Y ( t ) = ˜ Y cos(2 πf t ) + ˜ Y sin(2 πf t ) , (14)where we obtain the N × N matrices ˜ X , ˜ X , ˜ Y , and ˜ Y by substituting Eq. (14) into Eqs. (12) and (13) and then solving theresulting system of linear equations.We use root mean square (RMS) quantities as our principal diagnostic for evaluating our results, such as in our bifurcationdiagrams. Most commonly, we compute the RMS of the velocity of the y -component of the center particle (i.e., ˙ y , ). In thiscase, RMS = (cid:115) (cid:82) T ˙ y , ( t ) dtT , where T = 1 /f is the period of the excitation frequency. We show a plot of the RMS of the linear state (14) in Fig. 4(c) asa function of the excitation frequency for a fixed amplitude of a = 0 . mV and an excitation angle of φ = π/ . The loneresonant peak above the cutoff point is close to the estimated defect frequency f d ≈ . Hz.In Fig. 4(d), we show a frequency sweep in our experiment for a = 4 mV and φ = π/ . We show the theoretical values of thecutoff frequency f c ≈ . Hz and defect frequency f d ≈ . Hz that we found in Sec. III B as vertical solid and dashed lines,respectively. We observe that the experimental resonant peak is close to the theoretical value. To obtain a cleaner resonant peak,we use an excitation amplitude that is large enough to overcome the noise of the system. One such amplitude is a = 4 mV. Aswe will see in Sec. IV, an excitation amplitude of a = 4 mV is already in the nonlinear regime of the system. C. Numerical Methods for the Computation of Nonlinear Localized Modes and Their Stability
For the remainder of our paper, we focus on how the presence of nonlinearity affects the “defect-induced” linear localizedmodes of the system [see, e.g., Fig. 4(c)]. We refer to these solutions, which are localized in space and periodic in time, asnonlinear localized modes (NLMs). Unlike in the linear equations, we are unable to obtain analytical solutions for NLMs usingthe ansatz (14). Therefore, we turn to numerical computations. We compute time-periodic orbits of Eq. (3) with period T = 1 /f to high precision by finding roots of the map G = x ( T ) − x (0) , where x ( T ) is the solution of Eq. (3) at time T with initialcondition x (0) . Additionally, x ∈ R N is the vector that results from reshaping the matrix with elements x m,n , y m,n , ˙ x m,n ,and ˙ y m,n into row vectors and concatenating them into a single vector. We obtain roots of the map G using a Jacobian-freeNewton–Krylov method [40] with an initial guess of our linear state (14). We compute bifurcation diagrams using a pseudo-arclength continuation algorithm [41] with the excitation frequency f or amplitude a as our continuation parameter. Once weobtain the branches of a bifurcation diagram, we determine the linear stability of each solution x by solving the variationalequations V (cid:48) = DG · V with initial condition V (0) = I , where I denotes the identity matrix and DG is the N × N Jacobian matrix of the right-hand side of Eq. (3) evaluated at the solution x [42]. We calculate the Floquet multipliers, which (a) (b) -1 -0.5 0 0.5 1 Re( ) -1-0.500.51 I m () (c) Time [s]
Frequency [Hz] (d)
FIG. 5: Nonlinear localized mode of Eq. (3) that we obtain using a Newton–Krylov method for a = 4 mV, φ = π/ , and f = 9 . Hz. (a)
A surface plot of the NLM. We show the RMS velocity of each magnet in the lattice. (b)
Intensity plot of the NLM, where the color intensitycorresponds to the RMS. (c)
Floquet multipliers σ (blue marks) that are associated with the NLM in the complex plane. We show the unitcircle in gray. All multipliers lie within (or on) the unit circle, indicating that this solution is stable. (d) In the top panel, we show local maximaof the time series of ˙ y , when evolving zero initial data (i.e., the initial values of all variables are equal to ) for a = 4 mV, φ = π/ , and f = 9 . Hz. We approach the value ˙ y , ≈ . mm/s of the stable NLM. (See the black line.) In the bottom panel, we show the Fouriertransform of the final seconds of the time series of ˙ y , . This reveals a single large peak at frequency f ≈ . Hz. are hereafter denoted by σ , for a solution by computing the eigenvalues of the matrix V ( T ) . If all of the Floquet multipliers ofa solution have an absolute value that is less than or equal to , we say that the solution is “linearly stable”. Otherwise, we saythat the solution is “unstable”. Clearly, the Floquet multipliers only give information about the spectral stability of the solutions,and marginal instabilities associated with unit Floquet multipliers and nonlinear instabilities are possible. Therefore, stability isverified through numerical simulations. In our bifurcation diagrams, solid blue segments correspond to stable parameter regionsand dashed red segments correspond to unstable regions. We compute the Floquet multipliers after we obtain a solution with theNewton–Krylov method to avoid repeatedly solving the large variational system. This computation would be necessary if wewere using a standard Newton method, because the Jacobian of the map G is V ( T ) − I . IV. MAIN RESULTSA. Numerical NLMs
Using the Newton–Krylov method that we described in Sec. III C, we obtain a time-periodic solution with f = 9 . Hz, a = 4 mV , and φ = π/ . Additionally, because f = 9 . > . ≈ f c , this solution is localized in space [see Fig. 5(a)]. The dominantpeak is at the center of the lattice and the next-largest amplitudes are for the magnets that are adjacent to the center magnet atangles π/ , π/ , π/ , and π/ [see Fig. 5(b)]. This is not surprising, because we are exciting the lattice along the φ = π/ direction. The Floquet multipliers that are associated with this solution each have a magnitude that is no more than , indicatingthat the solution is stable [see Fig. 5(c)]. Indeed, upon simulating Eq. (3) with initial values of all variables equal to (i.e.,“zero initial data”) and f = 9 . Hz, a = 4 mV , and φ = π/ , the dynamics approaches this stable NLM. [See the top panel ofFig. 5(d).] As expected, the Fourier transform of the corresponding time series is localized around a frequency f ≈ . Hz.The spatial decay of the tails of the NLM depends on which direction of observation one considers. For example, if onemeasures the RMS velocity of the magnets that lie along the θ = π/ direction, the decay appears to be exponential or faster.See the solid blue squares in Fig. 6(a), which shows the RMS velocity versus distance from the origin (following the θ = π/ direction) in a semilog plot for the NLM from Fig. 5 [i.e., for the NLM with f = 9 . Hz, a = 4 mV , and φ = π/ ]. This isconsistent with the spatial decay properties of breathers in continuous-space settings, such as the ones in the quintic Ginzburg–Landau equation that were studied in [43]. The tails of the breathers in [43] decay at rate e − br / √ r , where b > is a constant.The solid yellow circles exhibit a similar decay for the magnets along the θ = 0 direction for our solution above, although weobserve some modulation in the decay profile, in contrast to the θ = π/ case. Modulations in spatial decay have been studiedin other settings, such as in the biharmonic φ model [44]. We observe similar decay properties for an NLM with f = 9 . Hz, a = 5 . mV and φ = π/ [see Fig. 6(b)]. B. Experimental NLMs
In our experiments, it is difficult to initialize the system with predetermined positions and velocities. To obtain the NLM, weexcite the system with a small amplitude ( a = 1 mV) which we increase gradually to the value a = 4 mV over about minutes.Because we predict the NLM to be stable at the resulting parameter values, we record data after 90 periods of motion have Distance from origin [mm] -2 -1 R M S [ mm / s ] (a) Distance from origin [mm] -2 -1 R M S [ mm / s ] (b) Time [s] -10001000 5 10 15
Frequency [Hz] (c)
FIG. 6: (a)
Decay of the NLM in the θ = π/ (blue squares) and θ = 0 (yellow circles) directions for a drive amplitude of a = 4 mV, driveangle of φ = π/ , and drive frequency of f = 9 . Hz. We show the RMS velocity versus the distance to the origin of the lattice. We showexperimental results as open markers and the numerical results as solid markers that are connected by lines. (b)
The same as panel (a), butfor a drive amplitude of a = 5 . mV. (c) Time series (top panel) of the center particle of the experimental NLM and corresponding Fouriertransform (bottom panel) with φ = π/ , f = 9 . Hz, and a = 4 mV. Despite the presence of some noise, the solution is predominantlyperiodic in time. Indeed, the Fourier transform of the time signal is highly localized around the frequency f = 9 . Hz. elapsed once we attain the value a = 4 mV. We track the velocities at the center particle with a laser vibrometer (see Sec. II),and we record the time series of the magnet velocities using an oscilloscope. [See the top panel of Fig. 6(c).] As expected, weobtain dynamics that are periodic in time, as we can see not only with the time series but also via its Fourier transform [whichwe show in the bottom panel of Fig. 6(c)]. We obtain similar experimental results for an amplitude of a = 5 . mV.To examine the spatial decay of the experimental NLM, we record the positions of the magnets in half of the lattice usinga digital camera (see Sec. II). By numerically differentiating the positions, we obtain an estimate for the velocities of thesemagnets. We were unable to do a complete full-field realization, because the DIC loses track of some magnets (if, e.g., themagnets begin to spin). However, we captured enough data to compute the decay along the two primary directions ( θ = π/ and θ = 0 ) that we examined in our numerical NLMs. The open blue squares of Fig. 6(a) show the decay along the θ = π/ direction and the open yellow circles show the decay along the θ = 0 direction for data that we obtained with f = 9 . Hz, a = 4 mV, and φ = π/ . The horizontal dashed line is our estimated mean value of the noise (see Sec. II). We show the experimentaldecay rates along with our numerical results. Although the numerical values overestimate the RMS velocity, the agreement isstill reasonable, especially for the center magnet. We find similar decay properties in our experiment with f = 9 . Hz, a = 5 . mV , and φ = π/ . [See the open markers of Fig. 6(b).]Recall that we do not tune the numerical results to fit the experimentally obtained NLM solution. Instead, we determine eachof the parameter values beforehand, as described in Sec. II. C. Frequency Continuation
In Figs. 4(d) and 6(a,b), we demonstrate that our model (3) agrees reasonably well with our experimental data. We nowconduct a series of numerical computations, with parameter continuation (see Sec. III C for a description of our procedure),to get a better sense of the role of nonlinearity in Eq. (3) and its interplay with the disorder (at the central magnet) and thediscreteness of the model. We return to our experiments in Sec. IV D to see what nonlinear effects we are able to capture in thelaboratory.We first perform continuation with respect to the excitation frequency f for a fixed excitation angle φ = π/ for various valuesof the excitation amplitude a . We thereby generate nonlinear analogs of the linear resonant peak that we showed in Fig. 4(c). InFig. 7(a), we show frequency continuations for our two drive amplitudes, a = 4 mV and a = 5 . , of the NLMs from Fig. 5 andFigs. 6(a,b). From comparing these frequency continuations to the linear case in Fig. 4(c), we see that the nonlinearity deformsthe peak, which becomes narrower and starts to bend towards higher frequencies. The nonlinearity also destabilizes the solutionsat some critical frequency; this occurs at f ≈ . Hz for a = 4 mV and at f ≈ . Hz for a = 5 . mV. We therefore seethat the NLM in Fig. 6(c) is unstable. Although our numerical computations predicted this NLM to be unstable, it is notable thatwe are able to access it in our experiments. Although this seems to imply that our theory is inconsistent with our experimentsfor the parameter values f ≈ . Hz and a = 5 . mV, the instability of the NLM for these parameter values is rather weak(with max i ( | σ i | ) ≈ . ). We observe instability after many periods when we perturb the NLM along the eigenvector that isassociated to the unstable Floquet multiplier σ ≈ . . See the top panel of Figure 7(b). However, if we initialize the dynamicswith zero initial data, we approach and stay close to the NLM solution that we obtained via a Newton–Krylov method even after200 periods of motion. See the bottom panel of Figure 7(b). This suggests that solutions with weak instabilities can still attractnearby points in phase space, at least initially, and that our numerical prediction for f ≈ . Hz and a = 5 . mV is consistent0 Frequency [Hz] R M S [ mm / s ] a = . m V a = m V (a) Time [s]
Time [s] (b)
Frequency [Hz] R M S [ mm / s ] a = . m V a = m V (c) FIG. 7: Frequency continuation of NLMs. (a)
Continuation with respect to frequency for fixed excitation angle φ = π/ for the twoamplitudes, a = 4 mV and a = 5 . mV, that we considered in Fig. 6. (b) In the top panel, we plot the local maxima of the time seriesof ˙ y , when we perturb the NLM with f = 9 . Hz and a = 5 . mV along the eigenvector that is associated with the largest-magnitudeFloquet multiplier. The size of this perturbation is equal to 5% of the amplitude of the solution. We show the corresponding value of the localmaximum of the NLM that we obtain via a Newton–Krylov method as the solid black line. In the bottom panel, we plot local maxima of thetime series when evolving zero initial data with a fixed frequency f = 9 . Hz and amplitude drive a = 5 . . (c) The same as panel (a), but withcolor intensity corresponding to the magnitude of the largest Floquet multiplier. Except near the peaks of the resonant curves, the instabilitiesare fairly weak.
Frequency [Hz] R M S [ mm / s ] a = m V a = m V a = m V a = m V (a) Frequency [Hz] R M S [ mm / s ] a = m V a = m V a = m V a = m V (b) Frequency [Hz] R M S [ mm / s ] a = m V a = m V a = . m V a = m V (c) FIG. 8: (a)
Frequency continuation with excitation angle φ = π/ for excitation amplitudes a = 1 mV, a = 2 mV, a = 4 mV, and a = 15 mV (b) Frequency continuation with φ = 0 for a = 1 mV, a = 2 mV, a = 4 mV, and a = 15 mV (c) Frequency continuation with φ = π/ for a = 1 mV, a = 1 . mV, a = 2 mV, and a = 3 mV. with our experimental observations. Figure 7(c) is the same as Figure 7(a), but now the color scale corresponds to the magnitudeof the maximum Floquet multiplier. This illustrates that the magnitude of the instability is weak throughout most of the solutionbranch. The instability is at its largest growth rate (with max i ( | σ i | ) ≈ . ) at the peak of the resonant curve.In Fig. 8(a), we show the gradual bending of the resonant peak for progressively larger excitation amplitudes with φ = π/ .In particular, for a = 15 mV, the peak bends so far that additional solutions at f ≈ . Hz emerge. However, these largeamplitude solutions are very unstable, and we were not able to access them either in our direct numerical simulations or in ourexperiments. Indeed, as we will discuss in Sec. IV D, we observe different types of dynamics at such excitation amplitudes. Wecan also tune the excitation angle φ and thereby deform the resonant peak in a different way. For example, when we fix theexcitation angle to φ = 0 (i.e., an excitation along n = 0 ), the resonant curves are qualitatively similar to those for φ = π/ for excitation amplitudes a = 1 , a = 2 , and a = 4 [see Fig. 8(b)], although the stability properties are slightly different. Forthe large excitation of a = 15 mV, the resonant curve bends even further toward higher frequencies. Even greater qualitativedifferences occur for φ = π/ (i.e., an excitation along m = 0 ); see Fig. 8(c). In this case, for small-amplitude excitations, theresonant peak has a unimodal shape, as expected. However, as we consider gradually larger excitation amplitudes, an additionalpeak begins to emerge from the main solution branch, leading to a two-humped profile in the dependence of the RMS velocity onthe frequency. This deformation is noticeable for relatively small excitation amplitudes. Specifically, we observe the existence ofmultiple solutions even for excitation amplitudes that are as small as a = 3 mV, where there is bifurcation at frequency f ≈ . Hz.The bifurcation diagram that we obtain when we use the excitation angle φ = π/ is more representative of “typical” anglesthan the special cases φ = π/ and φ = 0 . For example, even by decreasing the angle slightly from φ = π/ to φ = 89 π/ ,we observe the additional branch in the frequency continuation [see Fig. 9(a)]. We show a plot of an NLM at frequency f = 9 . Frequency [Hz] R M S [ mm / s ] = / = / (a) (b) (c) FIG. 9: (a)
Frequency continuation with an excitation amplitude of a = 3 mV for two excitation angles, φ = 89 π/ and φ = π/ . Ourcontinuations for φ = 89 π/ and φ = π/ are indistinguishable when we are outside the parameter region in which the φ = 89 π/ continuation has an additional branch of NLM solutions. (b) Profile of the NLM for f = 9 . Hz that belongs to the lower branch (i.e., thesmaller-amplitude NLM branch) of the φ = 89 π/ continuation. (b) Profile of the NLM for f = 9 . Hz that belongs to the upper branchof the φ = 89 π/ continuation. Drive amplitude [mV] R M S [ mm / s ] (a) -3 -3 (b) Drive amplitude [mV] R M S [ mm / s ] (c) FIG. 10: (a)
Our upsweep (dashed gray curve) and downsweep (solid black curve) of the drive amplitude a in our experiment. (b) Themagnitude of the Fourier transform of the time series for the velocity of the center particle normalized by the height of the peak at f ≈ . Hzfor (top panel) our experiment and (bottom panel) numerical computation for an excitation amplitude a = 5 . mV. The dashed gray curve isthe small-amplitude state (in the form of an NLM) that we obtain from the upsweep, and the solid black curve is the large-amplitude state thatwe obtain from the downsweep. (c) The upsweep (dashed gray curve) and downsweep (solid black curve) of the amplitude in our numericalsimulations.
Hz that belongs to the main branch (i.e., the branch with smaller-amplitude NLMs) of the φ = 89 π/ continuation in Fig. 9(b).It has a similar profile to the NLM in Fig. 5. We show a plot of an NLM from the additional (i.e., larger-amplitude) branch at thesame frequency f = 9 . Hz in Fig. 9(c). The solutions along this branch have secondary amplitudes in the − π/ direction. D. Drive-Amplitude Sweeps
We now return to the effect of large-amplitude excitations for the parameter set — namely, φ = π/ and f ≈ . Hz — inour laboratory experiments. Our bifurcation analysis revealed that the NLM solution at this parameter set destabilizes for largeramplitudes, although sometimes the instability is so weak that the NLMs are effectively stable on short enough time scales [seeFig. 7(b)]. We also observed that a large-amplitude branch of NLM solutions emerges at the parameter values φ = π/ and f ≈ . Hz for sufficiently large excitation amplitudes [see Fig. 8(a)].To study the dynamics at larger amplitudes in experiments, we initialize the system with a small-amplitude excitation (of a = 1 mV) and gradually increase the amplitude in increments of . mV. For each step, we run the system for 90 periods, whichallows sufficient time to settle to a steady state if there is one. We record the RMS of the velocity of the center magnet for the final5 seconds. We call this procedure an amplitude “upsweep”. We use an analogous procedure when we start with a large excitationamplitude, which we gradually decrease using the same increments. We call this procedure an amplitude “downsweep”. Weshow our experimental results for the upsweep and downsweep in Fig. 10(a). For sufficiently small amplitudes (specifically,for a (cid:47) . mV) both the upsweep and downsweep approach the same NLM, suggesting that there is a single stable branch2of NLMs for a (cid:47) . mV. However, for a (cid:39) . mV, there appear to be two different states; we obtain the small-amplitudestates when we perform an upsweep and the large-amplitude states when we perform a downsweep. The small-amplitude stateshave the form of an NLM. The experimental result in Fig. 6(b) is an example of the small-amplitude state for a = 5 . mV.The large-amplitude states are also localized, but they are no longer periodic in time. An inspection of the Fourier transform ofthe large-amplitude state for a = 5 . mV reveals other peaks in the spectrum (in addition to a peak at the excitation frequency f ≈ . Hz). In the top panel of Fig. 10(b), we show the Fourier transform of both the large-amplitude state and the small-amplitude state. The small-amplitude state (i.e., the NLM) has no peaks for lower frequencies, whereas the large-amplitude statehas peaks at approximately f = 8 . Hz and f = 6 . Hz; this is suggestive of quasiperiodic behavior.We obtain qualitatively similar results when we perform analogous upsweeps and downsweeps in numerical computations.We also observe the emergence of two states in these computations [see Fig. 10(c)]. In our computations, the large-amplitudestate departs from the branch of NLMs for excitation amplitudes that are slightly larger (specifically, for a (cid:39) . mV) than in ourexperiments. The amplitude a ≈ . mV is roughly where the numerical NLM branch destabilizes. As in our experiments, theselarge-amplitude states are not periodic in time. One can also observe the presence of secondary peaks in their Fourier transformsin our numerical solutions, although the locations of these peaks are slightly different than in our experiments. See the bottompanel of Fig. 10(b). Although these numerical large-amplitude states have features that are similar to those of time-quasiperiodicstates (given the multiple incommensurate peaks in the Fourier transform), it is also possible that these large-amplitude statesare weakly chaotic. One issue is that we were unable to detect asymptotically stable time-quasiperiodic orbits for parametervalues that correspond to the experiments. If such solutions were dynamic attractors, it would be straightforward to classify thelarge-amplitude states as time-quasiperiodic ones by plotting Poincar´e sections of the orbits.To clarify the nature of the large-amplitude states, we modify the parameter values slightly to obtain stable time-quasiperiodicsolutions. We perform amplitude upsweeps and downsweeps for the parameter set φ = π/ , f = 9 . Hz, and M b = 125 g. Weuse a different drive frequency from our prior calculations, because the smaller mass M b = 125 g leads to a cutoff frequencyof f c ≈ . Hz and a defect frequency of f d ≈ . Hz. The amplitude sweeps with these parameter values lead to a well-defined large-amplitude branch of solutions that bifurcates from the main branch of periodic ones (the NLMs) [see Fig. 11(a)].In Fig. 11(a), we show three solid markers to point out the locations of three solutions: a large-amplitude state that appears to beeither time-quasiperiodic or time-chaotic in black, the NLM (in gray), and a stable time-quasiperiodic orbit (in red). We showa plot of a projection of the Poincar´e section in the ( y, ˙ y ) plane in Fig. 11(b) for the two states that are not time-periodic. Theorbit in the bottom panel reveals a well-defined invariant curve, illustrating the quasiperiodic nature of the solution. We show theFourier transforms of these two non-periodic states in Fig. 11(c). Both have a secondary peak in the spectrum, demonstratingthat the solutions are indeed non-periodic in time (because of the incommensurate peaks in the frequency spectra). Laboratoryexperiments for this modified parameter set yield similar results. In particular, there is a well-defined large-amplitude branch ofsolutions that bifurcate from a branch of time-periodic solutions [see Fig. 11(d,e)]. The Fourier transform of one of the large-amplitude states also has a secondary peak in the spectrum, an indication that the state is nearly quasiperiodic. For simplicity,we henceforth call them simply “quasiperiodic”.Because the numerical computations and experiments with the modified mass M b = 125 g reveal the existence of time-quasiperiodic orbits that bifurcate from the main branch of time-periodic NLMs, it is reasonable to conclude that these quasiperi-odic solutions persist when we continue the parameters to the original parameter set M b = 138 . g and f d = 9 . . This suggeststhat the large-amplitude branch in Fig. 10(a,c) consists of time-quasiperiodic solutions. V. CONCLUSIONS
We have demonstrated, both experimentally and numerically, the existence of nonlinear localized modes in a 2D hexagonallattice of repelling magnets. By exploring the effects of nonlinearity numerically using frequency continuation and experimen-tally using amplitude sweeps, we revealed the emergence of both time-periodic NLMs and time-quasiperiodic localized states.We have also established that our experimental setup is a viable approach for fundamental studies in nonlinear lattice systemsthat go beyond what can occur in 1D chains. We found that the smaller-amplitude NLMs that we considered are stable, whereasprogressively larger excitation amplitudes led to instabilities and more complicated dynamics, including time-quasiperiodic andpotentially time-chaotic behavior. We also explored the anisotropy of the hexagonal lattice by considering different excitationangles and examining the nature and decay of the states along these angles.Our work paves the way for many future studies. For example, although our parameter continuation in frequency revealedseveral families of solutions, there are undoubtedly — given the complexity of the studied system — several other ones (includingpossibly exotic ones) to discover. Other avenues for future work include the study of refined models — such as ones that accountfor nonlinear damping (or, more generally, a more elaborate form of damping [45–47]), rotational effects (which are oftenconsidered to be important [48, 49]), and/or long-range interactions [37] — of our lattice system. Each of these aspects will addelements of complexity, but they also may lead to other types of interesting dynamics, such as the possibility of breather solutionswith algebraically decaying tails of waves in space. It is also certainly possible that the inclusion of rotational effects and/ormore sophisticated damping models may help improve matches with laboratory experiment. Such models have an associated3
Drive amplitude [mV] R M S [ mm / s ] (a) -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0-60-40-20-3.5 -3 -2.5 -2 -1.5 -1 -0.5 020304050 (b) Frequency (c)
Drive amplitude [mV] R M S [ mm / s ] (d) -5 0 5 Time (s) -1000100-5 0 5
Time (s) -1000100 (e)
Frequency [Hz]
Frequency [Hz] (f)
FIG. 11: (a)
The upsweep (dashed gray curve) and downsweep (solid black curve) for our numerical simulations with M b = 125 g and f = 9 . Hz. Three states exist when we use an excitation amplitude of a = 9 . mV. Two of the states (the black and red dots) are notperiodic in time, and the other (the gray dot) is time-periodic. We obtain the state that is indicated by the red dot by simulating zero initial datafor 90 periods of motion. (b) In the top panel, we show a projection of the Poincar´e section of the solution that is represented by the black dotin panel (a). In the bottom panel, we show a projection of the Poincar´e section of the solution that is represented by the red dot in panel (a). (c)
The magnitude of the Fourier transforms of the solutions in panel (b). (d)
The upsweep (dashed gray curve) and downsweep (solid blackcurve) for our experiments. (e)
The gray curve is the experimental time series of the small-amplitude state that we obtain from the upsweep foran excitation amplitude of a = 8 mV. The solid black curve is the time series of the large-amplitude state that we obtain from the downsweepfor an excitation amplitude of a = 8 mV. (f) The magnitude of the Fourier transform of the time series of panel (e) normalized by the heightof the peak at f = 9 . . The gray curve corresponds to the small-amplitude state, and the solid black curve corresponds to the large-amplitudestate. cost of being more complicated and hence more cumbersome to analyze and simulate. Our attempt in the present paper hasbeen to explore the principal features of the interplay of discreteness, local disorder, and nonlinearity in a hexagonal lattice ofmagnets. Breathers in heterogenous hexagonal magnetic lattices (e.g., ones with a repeating pattern of two masses) may lead tothe existence of intrinsic localized modes and are also worthy of future study. In that context, the study of band gaps, instabilities,and nonlinear modes and their propagation is another topic of substantial ongoing interest. Acknowledgements
The present paper is based on work that was supported by the US National Science Foundation under Grant Nos. DMS-1615037 (CC), DMS-1809074 (PGK), and EFRI-1741565 (CD). AJM acknowledges support from the Agencia Nacional deInvestigaci´on y Desarrollo de Chile (ANID) under Grant No. 3190906. EGC thanks Bowdoin College for the kind hospitalitywhere the initial stages of this work were carried out. PGK also acknowledges support from the Leverhulme Trust via a VisitingFellowship and thanks the Mathematical Institute of the University of Oxford for its hospitality during part of this work. We givespecial thanks to Bowdoin undergraduates Ariel Gonzales, Patrycja Pekala, Anam Shah, and Steven Xu for help with simulationsand data management.4
Appendix
To derive the external force that the wire exerts on a magnet, we first define our coordinate system. We choose a set oforthogonal unit vectors ˆ r , ˆ z , and ˆ s that are centered on the wire and oriented such that the wire is aligned with the ˆ s axis [seeFig. 12(a)]. We model the magnetic moment µ of the magnet using the Gilbert model of a magnetic dipole [50], F wire = ( µ · ∇ ) B : (15)which describes the force that acts on the dipole due to the magnetic field B . In our setup, the wire carries an electric current I that generates the magnetic field B ( r, z ) = B r ( r, z )ˆ r + B z ( r, z )ˆ z . Evaluating B at the position − h ˆ z + r ˆ r of the magnet yields B ( r, − h ) = Iµ π √ h + r (sin θ ˆ r + cos θ ˆ z ) , (16)where µ is the magnetic permeability and θ is the angle between B and µ . We are interested only in the ˆ r component of theforce, because we are assuming that we can neglect the dynamics along the ˆ z axis. Therefore, inserting Eq. (16) into Eq. (15)and taking µ = M ˆ z , we obtain F wire = M ∂ r B z ( r, z )ˆ r = Iµ M π h − r ( h + r ) ˆ r , (17)which corresponds to Eq. (1). In our coordinates, we place the wire in an orientation in the plane that is spanned by the latticebasis vectors e = (1 , and e = (1 / , √ / . However, it is straightforward to write Eq. (17) in coordinates in the { e , e } basis by including a parameter φ that accounts for the excitation angle. For instance, for the central magnet, we may write F ext0 , = Iµ M π cos( φ ) h − x , ( h + x , ) sin( φ ) h − y , ( h + y , ) , (18)which corresponds to Eq. (4). (a) (b) t [s] F o r c e M agn i t ude [ N ] -4 FIG. 12: (a)
Schematic of the interaction between a single magnet and a wire in the (ˆ r, ˆ z ) plane, which is orthogonal to the direction of thewire. We use “E.P.” to denote the equilibrium position of the magnet. (b) Squared magnitudes of the forces that result from lattice interactions(solid blue curve) and the external drive from the wire (dashed red curve) for the NLM from Fig. 5(a) during one period of motion. Specifically,we plot | F lattice0 , ( t ) | using the solid blue curve and | F ext0 , ( t ) | using the dashed red curve. Although the form of the force that describes the external drive is specific to our experimental setup, the dynamics aredominated by the lattice forces. For example, in Fig. 12(b), we compare the forces that result from the wire (see Eq. (18)) andthe force from the lattice for the NLM of Fig. 5(a). The lattice forces are F lattice m,n = − F ( q m +1 ,n − q m,n ) − F ( q m,n +1 − q m,n ) + F − ( q m,n − q m − ,n +1 )+ F ( q m,n − q m − ,n ) + F ( q m,n − q m,n − ) − F − ( q m +1 ,n − − q m,n ) (19)However, it is also important to acknowledge that the effective “defect” that is produced by the force (19) is responsible for the5presence of the corresponding linear defect frequency and hence for the associated NLMs in the presence of nonlinearity. [1] F. Lederer, G. I. Stegeman, D. N. Christodoulides, G. Assanto, M. Segev, and Y. Silberberg, “Discrete solitons in optics,” Phys. Rep. ,vol. 463, p. 1, 2008.[2] P. Binder, D. Abraimov, A. V. Ustinov, S. Flach, and Y. Zolotaryuk, “Observation of breathers in Josephson ladders,”
Phys. Rev. Lett. ,vol. 84, p. 745, 2000.[3] E. Tr´ıas, J. J. Mazo, and T. P. Orlando, “Discrete breathers in nonlinear lattices: Experimental detection in a Josephson array,”
Phys. Rev.Lett. , vol. 84, p. 741, 2000.[4] L. Q. English, M. Sato, and A. J. Sievers, “Modulational instability of nonlinear spin waves in easy-axis antiferromagnetic chains. II.influence of sample shape on intrinsic localized modes and dynamic spin defects,”
Phys. Rev. B , vol. 67, p. 024403, 2003.[5] U. T. Schwarz, L. Q. English, and A. J. Sievers, “Experimental generation and observation of intrinsic localized spin wave modes in anantiferromagnet,”
Phys. Rev. Lett. , vol. 83, p. 223, 1999.[6] B. I. Swanson, J. A. Brozik, S. P. Love, G. F. Strouse, A. P. Shreve, A. R. Bishop, W.-Z. Wang, and M. I. Salkola, “Observation ofintrinsically localized modes in a discrete low-dimensional material,”
Phys. Rev. Lett. , vol. 82, p. 3288, 1999.[7] M. Peyrard, “Nonlinear dynamics and statistical physics of DNA,”
Nonlinearity , vol. 17, p. R1, 2004.[8] J. Bajars, J. C. Eilbeck, and B. Leimkuhler,
Numerical Simulations of Nonlinear Modes in Mica: Past, Present and Future , vol. 221,p. 35. Cham, Switzerland: Springer International Publishing, 2015.[9] O. Morsch and M. Oberthaler, “Dynamics of Bose–Einstein condensates in optical lattices,”
Rev. Mod. Phys. , vol. 78, p. 179, 2006.[10] S. Flach and A. Gorbach, “Discrete breathers: Advances in theory and applications,”
Phys. Rep. , vol. 467, p. 1, 2008.[11] P. G. Kevrekidis, “Non-linear waves in lattices: Past, present, future,”
IMA J. Appl. Math. , vol. 76, pp. 389–423, 2011.[12] S. V. Dmitriev, E. A. Korznikova, Y. A. Baimova, and M. G. Velarde, “Discrete breathers in crystals,”
Physics-Uspekhi , vol. 59, pp. 446–461, may 2016.[13] E. Fermi, J. Pasta, and S. Ulam, “Studies of Nonlinear Problems. I.,” (Los Alamos National Laboratory, Los Alamos, NM, USA) , vol. Tech.Rep., pp. LA–1940, 1955.[14] G. Gallavotti,
The Fermi–Pasta–Ulam Problem: A Status Report . Heidelberg, Germany: Springer-Verlag, 2008.[15] V. F. Nesterenko,
Dynamics of Heterogeneous Materials . Heidelberg, Germany: Springer-Verlag, 2001.[16] C. Chong and P. G. Kevrekidis,
Coherent Structures in Granular Crystals: From Experiment and Modelling to Computation and Mathe-matical Analysis . Heidelberg, Germany: Springer-Verlag, 2018.[17] Y. Starosvetsky, K. R. Jayaprakash, M. A. Hasan, and A. F. Vakakis,
Topics On The Nonlinear Dynamics And Acoustics Of OrderedGranular Media . World Scientific, Singapore, 2017.[18] C. Chong, M. A. Porter, P. G. Kevrekidis, and C. Daraio, “Nonlinear coherent structures in granular crystals,”
J. Phys. Cond. Matt. ,vol. 29, p. 413003, 2017.[19] M. Moler´on, A. Leonard, and C. Daraio, “Solitary waves in a chain of repelling magnets,”
J. App. Phys. , vol. 115, no. 18, p. 184901,2014.[20] A. Mehrem, N. Jim´enez, L. J. Salmer´on-Contreras, X. Garc´ıa-Andr´es, L. M. Garc´ıa-Raffi, R. Pic´o, and V. J. S´anchez-Morcillo, “Nonlineardispersive waves in repulsive lattices,”
Phys. Rev. E , vol. 96, p. 012208, 2017.[21] M. Serra-Garcia, M. Moler´on, and C. Daraio, “Tunable, synchronized frequency down-conversion in magnetic lattices with defects,”
Phil. Trans. Roy. Soc. A: Math. Phys. Eng. Sci. , vol. 376, no. 2127, p. 20170137, 2018.[22] J. L. Mar´ın, J. C. Eilbeck, and F. M. Russell, , vol. 542, p. 293. Heidelberg, Germany: Springer-Verlag,2000.[23] L. Q. English, F. Palmero, J. F. Stormes, J. Cuevas, R. Carretero-Gonz´alez, and P. G. Kevrekidis, “Nonlinear localized modes in two-dimensional electrical lattices,”
Phys. Rev. E , vol. 88, p. 022912, 2013.[24] S. Vladimirov and K. Ostrikov, “Dynamic self-organization phenomena in complex ionized gas systems: New paradigms and technolog-ical aspects,”
Phys. Rep. , vol. 393, no. 3, pp. 175–380, 2004.[25] V. Koukouloyannis, P. G. Kevrekidis, K. J. H. Law, I. Kourakis, and D. J. Frantzeskakis, “Existence and stability of multisite breathers inhoneycomb and hexagonal lattices,”
J. Phys. A: Math. Theor. , vol. 43, no. 23, p. 235101, 2010.[26] J. A. D. Wattis,
Asymptotic Approximation of Discrete Breather Modes in Two-Dimensional Lattices , vol. 221, p. 179. Heidelberg,Germany: Springer-Verlag, 2015.[27] S. Flach, K. Kladko, and S. Takeno, “Acoustic breathers in two-dimensional lattices,”
Phys. Rev. Lett. , vol. 79, pp. 4838–4841, 1997.[28] J. L. Mar´ın, J. C. Eilbeck, and F. M. Russell, “Localized moving breathers in a 2D hexagonal lattice,”
Phys. Lett. A , vol. 248, no. 2,pp. 225–229, 1998.[29] B.-B. L¨u and T. Qiang, “Discrete gap breathers in a two-dimensional diatomic face-centered square lattice,”
Chinese Physics B , vol. 18,no. 10, p. 4393, 2009.[30] V. Koukouloyannis and I. Kourakis, “Discrete breathers in hexagonal dusty plasma lattices,”
Phys. Rev. E , vol. 80, p. 026402, 2009.[31] B.-F. Feng and T. Kawahara, “Discrete breathers in two-dimensional nonlinear lattices,”
Wave Motion , vol. 45, no. 1, pp. 68 – 82, 2007.[32] A. A. Maradudin, E. W. Montroll, and G. H. Weiss,
Theory of Lattice Dynamics in the Harmonic Approximation . New York City, NY,USA: Academic Press, 1963.[33] G. Theocharis, M. Kavousanakis, P. G. Kevrekidis, C. Daraio, M. A. Porter, and I. G. Kevrekidis, “Localized breathing modes in granularcrystals with defects,”
Phys. Rev. E , vol. 80, p. 066601, 2009.[34] N. Boechler, G. Theocharis, and C. Daraio, “Bifurcation based acoustic switching and rectification,”
Nat. Mater. , vol. 10, no. 9, p. 665, Euro. Phys. Lett. , vol. 101, p. 44003,2013.[36] C. Chong, A. Foehr, E. G. Charalampidis, P. G. Kevrekidis, and C. Daraio, “Breathers and other time-periodic solutions in an array ofcantilevers decorated with magnets,”
Math. Eng. , vol. 1, pp. 489–507, 2019.[37] M. Moler´on, C. Chong, A. J. Mart´ınez, M. A. Porter, P. G. Kevrekidis, and C. Daraio, “Nonlinear excitations in magnetic lattices withlong-range interactions,”
New J. Phys. , vol. 21, p. 063032, 2019.[38] S. Flach, “Breathers on lattices with long range interaction,”
Phys. Rev. E , vol. 58, pp. R4116–R4119, 1998.[39] R. F. Steidel,
An Introduction to Mechanical Vibrations . New York City, NY, USA: John Wiley & Sons, Inc., 1989.[40] C. T. Kelley,
Solving Nonlinear Equations with Newton’s Method . Philadelphia, PA, USA: Society for Industrial and Applied Mathemat-ics, 2003.[41] E. Doedel and L. S. Tuckerman,
Numerical Methods for Bifurcation Problems and Large-Scale Dynamical Systems . Heidelberg, Ger-many: Springer-Verlag, 2000.[42] M. W. Hirsch, S. Smale, and R. L. Devaney,
Differential Equations, Dynamical Systems, and an Introduction to Chaos . Amsterdam, TheNetherlands: Elsevier, 2004.[43] B. A. Malomed, “Potential of interaction between two- and three-dimensional solitons,”
Phys. Rev. E , vol. 58, pp. 7928–7933, Dec 1998.[44] R. J. Decker, A. Demirkaya, N. S. Manton, and P. G. Kevrekidis, “Kink–antikink interaction forces and bound states in a biharmonic φ model,” arXiv:2005.04523 , 2020.[45] A. Rosas, A. H. Romero, V. F. Nesterenko, and K. Lindenberg, “Observation of two-wave structure in strongly nonlinear dissipativegranular chains,” Phys. Rev. Lett. , vol. 98, p. 164301, 2007.[46] L. Vergara, “Model for dissipative highly nonlinear waves in dry granular systems,”
Phys. Rev. Lett. , vol. 104, p. 118001, 2010.[47] R. Carretero-Gonz´alez, D. Khatri, M. A. Porter, P. G. Kevrekidis, and C. Daraio, “Dissipative solitary waves in granular crystals,”
Phys.Rev. Lett. , vol. 102, p. 024102, 2009.[48] A. Merkel, V. Tournat, and V. Gusev, “Experimental evidence of rotational elastic waves in granular phononic crystals,”
Phys. Rev. Lett. ,vol. 107, p. 225502, 2011.[49] K. Vorotnikov, M. Kovaleva, and Y. Starosvetsky, “Emergence of non-stationary regimes in one- and two-dimensional models withinternal rotators,”
Phil. Trans. A: Math. Phys. Eng. Sci. , vol. 376, p. 20170134, 2018.[50] J. D. Jackson,