Nonlinear Localized Modes in Two-Dimensional Hexagonally-Packed Magnetic Lattices

We conduct an extensive study of nonlinear localized modes (NLMs), which are temporally periodic and spatially localized structures, in a two-dimensional array of repelling magnets. In our experiments, we arrange a lattice in a hexagonal configuration with a light-mass defect, and we harmonically drive the center of the chain with 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, we obtain good qualitative agreement between theory and experiments for a variety of dynamical behaviors. We find that the spatial decay is direction-dependent and that drive amplitudes along fundamental displacement axes lead to nonlinear resonant peaks in frequency continuations that are similar to those that occur in one-dimensional damped, driven lattices. However, driving along other directions leads to the creation of asymmetric NLMs that bifurcate from the main solution branch, which consists of symmetric NLMs. When we vary the drive amplitude, we observe such behavior both in our experiments and in our simulations. We also demonstrate that solutions that appear to be time-quasi-periodic bifurcate from the branch of symmetric time-periodic NLMs.

The 2D setting of the present work is a mechanical system in which each magnet has two displacement fields, which distinguishes it from many studies of scalar 2D lattices, such as those that describe electrical circuits [23]. Specifically, we examine a lattice of repelling magnets that are arranged in a hexagonal configuration. The choice of a hexagonal arrangement is motivated by the experimental setup, as hexagonal configurations are more robust structurally than other arrangements (such as square ones). At the center of the lattice is a light-mass defect, which introduces a localized defect mode into the spectrum of the linearization of the system. To excite the system experimentally, we drive the center of the lattice by a force that results from the current that flows along a wire that we suspend above the lattice. We model damping using a dashpot term. Putting everything together, 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 of breathers. Linear systems with an impurity or a defect (e.g., with a particle of lighter mass than the other particles) have isolated points in their spectra that lie above the spectral edge. These modes are called "defect localized modes" [32]. In the presence of 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 way for 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 nonlinear lattices. For example, such orbits have been identified in strongly nonlinear damped, driven granular chains [34,35], suggesting 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 we explore NLM profiles, spatial decay, parameter continuations, and nearly time-quasiperiodic orbits. We conclude and discuss future 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 attach a neodymium magnetic cylinder. We glue a thin piece of cover glass at the bottom to make the surface smoother and thereby improve the levitation of the particles. We build the defect particle, which is located in the center of the lattice, by directly attaching the magnet on the glass without the 3D-printed structure. This particle has a lighter mass and serves as a defect [see Fig. 1(b)]. The mean mass of a normal disc particle is 138.2 mg ± 3.1 mg (where we measure the standard deviation from a sample of 20 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 the particle at a height of 3 mm. We generate the AC current that flows through the wire from a Lock-In amplifier (SR860 500kHz DSP Lock-in Amplifier), and we amplify it with an audio amplifier (Topping TP22, class D). The equation that describes the force that the wire exerts on a magnet at distance r from it can be calculated as: where h is the height of the wire from the plane of floating discs, I is the wire current, µ 0 = 4π · 10 −7 N A −2 is the magnetic permeability, and M = 7.8 · 10 −3 Am 2 is the magnetic moment of the floating disc. See the appendix for the derivation of Eq. (1). We use harmonic excitations in our experiment, so the current through the wire is I(t) = aI 0 sin(2πf t), where f is the drive frequency (in Hz), a is the drive-voltage amplitude (in Volts), and I 0 = 0.1 A V −1 is the current per unit voltage that we measure in the wire. The magnets repel each other. In the ideal situation of a perfect dipole-dipole interaction, the magnetic force between two repelling magnets is where r is the distance (in meters) between the two center points of the magnets, p = −4, and A = 3µ 0 M 2 /2π. Although Eq. (2) is reasonable for large separation distances, we obtain better agreement by empirically determining A and p. Because the force between two magnetic dipoles is too small to measure directly, we create a pair of plastic plates, with 25 magnets attached to each plate. We position the plates to align each pair of cylindrical magnets from opposite plates through their radial directions. 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 that we can neglect interactions between magnets that are not aligned. The distance between a magnet on the first plate and the non-aligned magnets of the other plate is larger than 25 mm. As one can see in Fig. 1(c) the interaction force already tends to 0 for distances smaller that are significantly smaller than 25 mm. Consequently, the measured force is approximately equal to the sum of the repulsive force of the 25 isolated magnet pairs. We fit the data using Eq. (2), which yields p ≈ −4.2 and A ≈ 3.8 · 10 −12 N/m p [see Fig. 1 We monitor the motion of the central particle using a laser vibrometer (Polytec CLV-2534), and we record the remainder of the lattice using a digital camera (Point Grey GS3-U3-41C6C-C) with a frame rate of 90 fps. We analyze the images using digital-image-correlation (DIC) software (VIC-2D) to determine each particle's velocity. We inspect half of the lattice, as the cables that are connected to the driving wire block most of the system's other half [see Fig. 1(a)]. Due to imperfections at the bottom 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 particles start to rotate when they are levitated by the air that flows out of the air-bearing table. The image correlation software then loses track of them. We ignore these rotating discs in our subsequent analysis. To estimate the value of the damping coefficient, we excite the center particle and record the resulting temporal amplitude decay once we switch off the excitation. We fit the decay to an exponential function, which we then use to calculate the damping coefficient γ ≈ 10.52 · 10 −3 N s/m. The lattice particles are always in motion with at least small speeds, even in the absence of excitation. This is due to interactions with the air flow from the table and to imperfections (e.g., nonaxisymmetric mass distributions) of the particles. We use this motion to estimate the noise in the system. To evaluate the amount of noise, we record the lattice motion without excitation as a comparison. We summarize the values of the parameters in Table I.

A. 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 are essential for NLMs and also yield reasonable agreement with experimental data. It is in this spirit that we develop our model equations. After doing so (at the end of this section), we briefly discuss model simplifications.
One can express the distance between the magnet with index (m, n) and one of its nearest neighbors in terms of the displacements x m,n and y m,n of the magnets from their respective equilibrium positions. Once we determine the distance, we compute the resulting force using Eq. (2). Summing the forces from each of the six nearest neighbors and applying Newton's second law leads to the following equations of motion: When in equilibrium, the center-to-center distance between any two adjacent magnets is δ. The outer hexagonal boundary is the solid gray hexagon that encloses the lattice. On the boundary, the solid points represent fixed (i.e., immovable) magnets. There are w + 1 such magnets along each edge of the boundary. The solid circle represents the defect particle, which has index (m, n) = (0, 0).

The vector functions
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 account for damping. Using such a term has yielded reasonable agreement with experiments in other, similar lattices [19,36,37]. The quantity 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 a wire that is directly above the center of the lattice. The magnitude of the excitation is given by Eq. (1). Therefore, where φ is the angle of the excitation and F ext m,n = 0 when m = 0 and n = 0. In our experiments and in most of our numerical computations, the excitation angle is φ = π/2, so we excite only the y-component of the center magnet. We will also explore some other excitation angles. As we discuss in the appendix, the lattice forces dominate the dynamics. The wire has only a small effect on magnets that are beyond the center of the lattice. For example, at equilibrium, the force that is exerted on the center magnet 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 longrange effects can alter the structure of localized modes. For example, it was shown in [38] that the spatial decay of breathers can transition from exponential spatial decay to algebraic decay in lattices with algebraically decaying interaction forces (as is the case in our model) for lattices with sufficiently many sites. More recently, Molerón et al. [37] studied NLMs in a 1D magnetic lattice using a model with long-range interactions. Although the differences between long-range and nearest-neighbor lattices that were considered in [37] are detectable, they are still small. For example, at equilibrium, the force that is exerted on the center magnet by its nearest neighbors is one order of magnitude larger than that exerted by its next-nearest neighbor. (Compare Eq. (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. This leaves air resistance as the primary source of damping. Given the size of the magnets and velocities that we consider, we employ a linear dashpot [39]. We also assume that the magnets stay in a plane. We validate the many assumptions that we have made in formulating 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 the excitation amplitude a, frequency f , and angle φ. We will specify these in our various examples. In all cases, we examine a 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 beforehand using the procedures that we detailed in Sec. II.

B. Linear Analysis
We start with the basic linear theory of localized modes for our hexagonal magnetic lattice. We are particularly interested in modes with frequencies that lie above the cutoff frequency of the pass band. We first derive an analytical expression for the cutoff frequency, which is straightforward for an infinite-dimensional Hamiltonian system (i.e., with all integers m and n, along with a = 0 and γ = 0). We then numerically estimate the frequency of a linear mode that is associated with the defect for in the finite-dimensional Hamiltonian system. Finally, we compute linear localized modes in the associated finite-dimensional damped, driven system.
Assuming small strains, such that one can Taylor expand, where DF j is the Jacobian matrix of F j . Using this notation, we write the linearized equations of motion: where the entries of the Jacobian matrices are denoted by, whereδ ≡ Aδ p−1 . For a monoatomic system (in which all magnets are identical, such that M m,n = M b ), the linear system has plane-wave solutions where the wavenumbers k, and angular frequency ω = ω(k, ) satisfy the dispersion relationship where Fig. 3(a,b), we show contour plots of the two dispersion surfaces from Eq. (7). The dispersion curves along the edge of the irreducible Brillouin zone are shown in Fig. 3(c). The cutoff value of the pass band has the wavenumber pair (k, ) = (0, 2π/3), which is where the dispersion curve attains its maximum value. For the parameter values in Table I, the cutoff frequency is f c = ω(0, 2π/3) + /(2π) ≈ 8.77 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 a frequency above the cutoff frequency of the linear monoatomic system. With the light-mass defect at the center of the lattice, we write where (0, 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., the non-defect mass), and M δ < M b . We numerically compute the linear modes of the system with a mass defect, and we are thereby able to consider finite lattices. We use a hexagonal boundary with edge length w = 6 magnets [see Fig. 2(b)]. One can embed this lattice into a square matrix of size N × N , where N = 2w − 1 is the number magnets along the n = 0 line of the lattice. 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)th entry is y m,n (t). We enforce the fixed hexagonal boundaries by setting the displacements of magnets with indices (m, n) such that |m + n| > w to 0. We define the N × N matrix operators L α through where α ∈ {a, b, c, d}; the N × N tridiagonal matrix D has 1 entries on the super-diagonals and sub-diagonals, 2 entries along the diagonal, and 0 entries everywhere else; E is an N × N matrix with 1 entries along the super-diagonal and 0 entries everywhere else; and E T is the transpose of E. With these definitions, Eq. (6) becomes where M is an N × N matrix in which all entries except the (0, 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) =Ỹ e iωt , whereX andỸ are N × N time-independent matrices and One can cast Eq. (11) as a standard eigenvalue problem by letting λ = −ω 2 and unwrapping theX andỸ matrices into equivalent row vectors and reshaping the block matrix (with entries given by L α ) into a corresponding 2N 2 × 2N 2 matrix. One can then numerically solve the resulting eigenvalue problem to obtain 2N 2 eigenvalues and their corresponding modes. Using the 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 ≈ 9.17) lie above the cutoff frequency f c ≈ 8.77 Hz. We plot the value of f d as the horizontal dashed line in Fig. 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 approximate system: where A is an N × N matrix that has all 0 entries except for the single nonzero entry 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 where we obtain the N × N matricesX 1 ,X 2 ,Ỹ 1 , andỸ 2 by substituting Eq. (14) into Eqs. (12) and (13) and then solving the resulting system of linear equations. We use root mean square (RMS) quantities as our principal diagnostic for evaluating our results, such as in our bifurcation diagrams. Most commonly, we compute the RMS of the velocity of the y-component of the center particle (i.e.,ẏ 0,0 ). In this case, 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) as a function of the excitation frequency for a fixed amplitude of a = 0.01 mV and an excitation angle of φ = π/2. The lone resonant peak above the cutoff point is close to the estimated defect frequency f d ≈ 9.17 Hz. In Fig. 4(d), we show a frequency sweep in our experiment for a = 4 mV and φ = π/2. We show the theoretical values of the cutoff frequency f c ≈ 8.77 Hz and defect frequency f d ≈ 9.17 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. As we 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 localized modes of the system [see, e.g., Fig. 4(c)]. We refer to these solutions, which are localized in space and periodic in time, as nonlinear localized modes (NLMs). Unlike in the linear equations, we are unable to obtain analytical solutions for NLMs using the ansatz (14). Therefore, we turn to numerical computations. We compute time-periodic orbits of Eq. is the vector that results from reshaping the matrix with elements x m,n , y m,n ,ẋ m,n , andẏ m,n into row vectors and concatenating them into a single vector. We obtain roots of the map G using a Jacobian-free Newton-Krylov method [40] with an initial guess of our linear state (14). We compute bifurcation diagrams using a pseudoarclength continuation algorithm [41] with the excitation frequency f or amplitude a as our continuation parameter. Once we obtain the branches of a bifurcation diagram, we determine the linear stability of each solution x by solving the variational equations V = DG · V with initial condition V (0) = I, where I denotes the identity matrix and DG is the 4N 2 × 4N 2 Jacobian matrix of the right-hand side of Eq. (3) evaluated at the solution x [42]. We calculate the Floquet multipliers, which are hereafter denoted by σ, for a solution by computing the eigenvalues of the matrix V (T ). If all of the Floquet multipliers of a solution have an absolute value that is less than or equal to 1, we say that the solution is "linearly stable". Otherwise, we say that 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 is verified through numerical simulations. In our bifurcation diagrams, solid blue segments correspond to stable parameter regions and dashed red segments correspond to unstable regions. We compute the Floquet multipliers after we obtain a solution with the Newton-Krylov method to avoid repeatedly solving the large variational system. This computation would be necessary if we were using a standard Newton method, because the Jacobian of the map G is V (T ) − I.

A. Numerical NLMs
Using the Newton-Krylov method that we described in Sec. III C, we obtain a time-periodic solution with f = 9.3 Hz, a = 4 mV , and φ = π/2. Additionally, because f = 9.3 > 8.77 ≈ f c , this solution is localized in space [see Fig. 5(a)]. The dominant peak is at the center of the lattice and the next-largest amplitudes are for the magnets that are adjacent to the center magnet at angles π/3, 2π/3, 4π/3, and 5π/3 [see Fig. 5(b)]. This is not surprising, because we are exciting the lattice along the φ = π/2 direction. The Floquet multipliers that are associated with this solution each have a magnitude that is no more than 1, indicating that the solution is stable [see Fig. 5(c)]. Indeed, upon simulating Eq. (3) with initial values of all variables equal to 0 (i.e., "zero initial data") and f = 9.3 Hz, a = 4 mV , and φ = π/2, the dynamics approaches this stable NLM. [See the top panel of Fig. 5(d).] As expected, the Fourier transform of the corresponding time series is localized around a frequency f ≈ 9.3 Hz.
The spatial decay of the tails of the NLM depends on which direction of observation one considers. For example, if one measures the RMS velocity of the magnets that lie along the θ = π/3 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 θ = π/3 direction) in a semilog plot for the NLM from Fig. 5 [i.e., for the NLM with f = 9.3 Hz, a = 4 mV , and φ = π/2]. This is consistent 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 > 0 is a constant.
The solid yellow circles exhibit a similar decay for the magnets along the θ = 0 direction for our solution above, although we observe some modulation in the decay profile, in contrast to the θ = π/3 case. Modulations in spatial decay have been studied in other settings, such as in the biharmonic φ 4 model [44]. We observe similar decay properties for an NLM with f = 9.3 Hz, a = 5.5 mV and φ = π/2 [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, we excite the system with a small amplitude (a = 1 mV) which we increase gradually to the value a = 4 mV over about 3 minutes. Because we predict the NLM to be stable at the resulting parameter values, we record data after 90 periods of motion have To examine the spatial decay of the experimental NLM, we record the positions of the magnets in half of the lattice using a digital camera (see Sec. II). By numerically differentiating the positions, we obtain an estimate for the velocities of these magnets. We were unable to do a complete full-field realization, because the DIC loses track of some magnets (if, e.g., the magnets begin to spin). However, we captured enough data to compute the decay along the two primary directions (θ = π/3 and θ = 0) that we examined in our numerical NLMs. The open blue squares of Fig. 6(a) show the decay along the θ = π/3 direction and the open yellow circles show the decay along the θ = 0 direction for data that we obtained with f = 9.3 Hz, a = 4 mV, and φ = π/2. The horizontal dashed line is our estimated mean value of the noise (see Sec. II). We show the experimental decay rates along with our numerical results. Although the numerical values overestimate the RMS velocity, the agreement is still reasonable, especially for the center magnet. We find similar decay properties in our experiment with f = 9.3 Hz, a = 5.5 mV , and φ = π/2. [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 each of 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 now conduct 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 the discreteness of the model. We return to our experiments in Sec. IV D to see what nonlinear effects we are able to capture in the laboratory.
We first perform continuation with respect to the excitation frequency f for a fixed excitation angle φ = π/2 for various values of the excitation amplitude a. We thereby generate nonlinear analogs of the linear resonant peak that we showed in Fig. 4(c). In Fig. 7(a), we show frequency continuations for our two drive amplitudes, a = 4 mV and a = 5.5, of the NLMs from Fig. 5 and Figs. 6(a,b). From comparing these frequency continuations to the linear case in Fig. 4(c), we see that the nonlinearity deforms the peak, which becomes narrower and starts to bend towards higher frequencies. The nonlinearity also destabilizes the solutions at some critical frequency; this occurs at f ≈ 9.29 Hz for a = 4 mV and at f ≈ 9.31 Hz for a = 5.5 mV. We therefore see that the NLM in Fig. 6(c) is unstable. Although our numerical computations predicted this NLM to be unstable, it is notable that we are able to access it in our experiments. Although this seems to imply that our theory is inconsistent with our experiments for the parameter values f ≈ 9.31 Hz and a = 5.5 mV, the instability of the NLM for these parameter values is rather weak (with max i (|σ i |) ≈ 1.007). We observe instability after many periods when we perturb the NLM along the eigenvector that is associated to the unstable Floquet multiplier σ ≈ 1.007. See the top panel of Figure 7(b). However, if we initialize the dynamics with zero initial data, we approach and stay close to the NLM solution that we obtained via a Newton-Krylov method even after 200 periods of motion. See the bottom panel of Figure 7(b). This suggests that solutions with weak instabilities can still attract nearby points in phase space, at least initially, and that our numerical prediction for f ≈ 9.31 Hz and a = 5.5 mV is consistent  with our experimental observations. Figure 7(c) is the same as Figure 7(a), but now the color scale corresponds to the magnitude of the maximum Floquet multiplier. This illustrates that the magnitude of the instability is weak throughout most of the solution branch. The instability is at its largest growth rate (with max i (|σ i |) ≈ 2.15) 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 φ = π/2. In particular, for a = 15 mV, the peak bends so far that additional solutions at f ≈ 9.3 Hz emerge. However, these large amplitude solutions are very unstable, and we were not able to access them either in our direct numerical simulations or in our experiments. Indeed, as we will discuss in Sec. IV D, we observe different types of dynamics at such excitation amplitudes. We can also tune the excitation angle φ and thereby deform the resonant peak in a different way. For example, when we fix the excitation angle to φ = 0 (i.e., an excitation along n = 0), the resonant curves are qualitatively similar to those for φ = π/2 for excitation amplitudes a = 1, a = 2, and a = 4 [see Fig. 8(b)], although the stability properties are slightly different. For the large excitation of a = 15 mV, the resonant curve bends even further toward higher frequencies. Even greater qualitative differences occur for φ = π/3 (i.e., an excitation along m = 0); see Fig. 8(c). In this case, for small-amplitude excitations, the resonant peak has a unimodal shape, as expected. However, as we consider gradually larger excitation amplitudes, an additional peak begins to emerge from the main solution branch, leading to a two-humped profile in the dependence of the RMS velocity on the frequency. This deformation is noticeable for relatively small excitation amplitudes. Specifically, we observe the existence of multiple solutions even for excitation amplitudes that are as small as a = 3 mV, where there is bifurcation at frequency f ≈ 9.28 Hz.
The bifurcation diagram that we obtain when we use the excitation angle φ = π/3 is more representative of "typical" angles than the special cases φ = π/2 and φ = 0. For example, even by decreasing the angle slightly from φ = π/2 to φ = 89π/180, we observe the additional branch in the frequency continuation [see Fig. 9(a)]. We show a plot of an NLM at frequency f = 9.2  Hz that belongs to the main branch (i.e., the branch with smaller-amplitude NLMs) of the φ = 89π/180 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 the same frequency f = 9.2 Hz in Fig. 9(c). The solutions along this branch have secondary amplitudes in the −π/3 direction.

D. Drive-Amplitude Sweeps
We now return to the effect of large-amplitude excitations for the parameter set -namely, φ = π/2 and f ≈ 9.3 Hz -in our laboratory experiments. Our bifurcation analysis revealed that the NLM solution at this parameter set destabilizes for larger amplitudes, although sometimes the instability is so weak that the NLMs are effectively stable on short enough time scales [see Fig. 7(b)]. We also observed that a large-amplitude branch of NLM solutions emerges at the parameter values φ = π/2 and f ≈ 9.3 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 0.05 mV. For each step, we run the system for 90 periods, which allows 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 final 5 seconds. We call this procedure an amplitude "upsweep". We use an analogous procedure when we start with a large excitation amplitude, which we gradually decrease using the same increments. We call this procedure an amplitude "downsweep". We show our experimental results for the upsweep and downsweep in Fig. 10(a). For sufficiently small amplitudes (specifically, for a 4.5 mV) both the upsweep and downsweep approach the same NLM, suggesting that there is a single stable branch of NLMs for a 4.5 mV. However, for a 4.5 mV, there appear to be two different states; we obtain the small-amplitude states when we perform an upsweep and the large-amplitude states when we perform a downsweep. The small-amplitude states have the form of an NLM. The experimental result in Fig. 6(b) is an example of the small-amplitude state for a = 5.5 mV. The large-amplitude states are also localized, but they are no longer periodic in time. An inspection of the Fourier transform of the large-amplitude state for a = 5.5 mV reveals other peaks in the spectrum (in addition to a peak at the excitation frequency f ≈ 9.3 Hz). In the top panel of Fig. 10(b), we show the Fourier transform of both the large-amplitude state and the smallamplitude state. The small-amplitude state (i.e., the NLM) has no peaks for lower frequencies, whereas the large-amplitude state has peaks at approximately f = 8.9 Hz and f = 6.2 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-amplitude state departs from the branch of NLMs for excitation amplitudes that are slightly larger (specifically, for a 5.1 mV) than in our experiments. The amplitude a ≈ 5.1 mV is roughly where the numerical NLM branch destabilizes. As in our experiments, these large-amplitude states are not periodic in time. One can also observe the presence of secondary peaks in their Fourier transforms in our numerical solutions, although the locations of these peaks are slightly different than in our experiments. See the bottom panel of Fig. 10(b). Although these numerical large-amplitude states have features that are similar to those of time-quasiperiodic states (given the multiple incommensurate peaks in the Fourier transform), it is also possible that these large-amplitude states are weakly chaotic. One issue is that we were unable to detect asymptotically stable time-quasiperiodic orbits for parameter values that correspond to the experiments. If such solutions were dynamic attractors, it would be straightforward to classify the large-amplitude states as time-quasiperiodic ones by plotting Poincaré sections of the orbits.
To clarify the nature of the large-amplitude states, we modify the parameter values slightly to obtain stable time-quasiperiodic solutions. We perform amplitude upsweeps and downsweeps for the parameter set φ = π/2, f = 9.65 Hz, and M b = 125 g. We use a different drive frequency from our prior calculations, because the smaller mass M b = 125 g leads to a cutoff frequency of f c ≈ 9.01 Hz and a defect frequency of f d ≈ 9.36 Hz. The amplitude sweeps with these parameter values lead to a welldefined 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 be either time-quasiperiodic or time-chaotic in black, the NLM (in gray), and a stable time-quasiperiodic orbit (in red). We show a plot of a projection of the Poincaré section in the (y,ẏ) plane in Fig. 11(b) for the two states that are not time-periodic. The orbit in the bottom panel reveals a well-defined invariant curve, illustrating the quasiperiodic nature of the solution. We show the Fourier transforms of these two non-periodic states in Fig. 11(c). Both have a secondary peak in the spectrum, demonstrating that the solutions are indeed non-periodic in time (because of the incommensurate peaks in the frequency spectra). Laboratory experiments for this modified parameter set yield similar results. In particular, there is a well-defined large-amplitude branch of solutions that bifurcate from a branch of time-periodic solutions [see Fig. 11(d,e)]. The Fourier transform of one of the largeamplitude 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 timequasiperiodic orbits that bifurcate from the main branch of time-periodic NLMs, it is reasonable to conclude that these quasiperiodic solutions persist when we continue the parameters to the original parameter set M b = 138.2 g and f d = 9.3. This suggests that 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 hexagonal lattice of repelling magnets. By exploring the effects of nonlinearity numerically using frequency continuation and experimentally 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 systems that go beyond what can occur in 1D chains. We found that the smaller-amplitude NLMs that we considered are stable, whereas progressively larger excitation amplitudes led to instabilities and more complicated dynamics, including time-quasiperiodic and potentially time-chaotic behavior. We also explored the anisotropy of the hexagonal lattice by considering different excitation angles 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 revealed several families of solutions, there are undoubtedly -given the complexity of the studied system -several other ones (including possibly exotic ones) to discover. Other avenues for future work include the study of refined models -such as ones that account for nonlinear damping (or, more generally, a more elaborate form of damping [45][46][47]), rotational effects (which are often considered to be important [48,49]), and/or long-range interactions [37] -of our lattice system. Each of these aspects will add elements of complexity, but they also may lead to other types of interesting dynamics, such as the possibility of breather solutions with algebraically decaying tails of waves in space. It is also certainly possible that the inclusion of rotational effects and/or more sophisticated damping models may help improve matches with laboratory experiment. Such models have an associated cost of being more complicated and hence more cumbersome to analyze and simulate. Our attempt in the present paper has been to explore the principal features of the interplay of discreteness, local disorder, and nonlinearity in a hexagonal lattice of magnets. Breathers in heterogenous hexagonal magnetic lattices (e.g., ones with a repeating pattern of two masses) may lead to the 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.
To derive the external force that the wire exerts on a magnet, we first define our coordinate system. We choose a set of orthogonal unit vectorsr,ẑ, andŝ that are centered on the wire and oriented such that the wire is aligned with theŝ axis [see Fig. 12(a)]. We model the magnetic moment µ of the magnet using the Gilbert model of a magnetic dipole [50], 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)ẑ. Evaluating B at the position −hẑ + rr of the magnet yields B(r, −h) = Iµ 0 2π √ h 2 + r 2 (sin θr + cos θẑ) , where µ 0 is the magnetic permeability and θ is the angle between B and µ. We are interested only in ther component of the force, because we are assuming that we can neglect the dynamics along theẑ axis. Therefore, inserting Eq. (16) into Eq. (15) and taking µ = Mẑ, we obtain which corresponds to Eq. (1). In our coordinates, we place the wire in an orientation in the plane that is spanned by the lattice basis vectors e 1 = (1, 0) and e 2 = (1/2, √ 3/2). However, it is straightforward to write Eq. (17) in coordinates in the {e 1 , e 2 } basis by including a parameter φ that accounts for the excitation angle. For instance, for the central magnet, we may write which corresponds to Eq. (4).  12: (a) Schematic of the interaction between a single magnet and a wire in the (r,ẑ) plane, which is orthogonal to the direction of the wire. 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 lattice 0,0 (t)| 2 using the solid blue curve and |F ext 0,0 (t)| 2 using the dashed red curve.
Although the form of the force that describes the external drive is specific to our experimental setup, the dynamics are dominated by the lattice forces. For example, in Fig. 12(b), we compare the forces that result from the wire (see Eq. (18)) and the force from the lattice for the NLM of Fig. 5(a). The lattice forces are F lattice m,n = −F 0 (q m+1,n − q m,n ) − F 1 (q m,n+1 − q m,n ) + F −1 (q m,n − q m−1,n+1 ) +F 0 (q m,n − q m−1,n ) + F 1 (q m,n − q m,n−1 ) − F −1 (q m+1,n−1 − q m,n ) However, it is also important to acknowledge that the effective "defect" that is produced by the force (19) is responsible for the