Nonlinearities in dynamical systems give rise to a variety of interesting physical eects, many of which have the potential to be harnessed for technological applications.
One of the most well-known examples is the generation of an optical frequency comb [1], in which a coherent and short optical laser pulse train results in a discrete spectrum comprising millions of optical modes. In the context of communications, compact and phase-synchronous frequency combs can be produced in microresonators taking advantage of Kerr nonlinearity [2] or Pockels eect [3] which self-modulate an electromagnetic wave. More recently, frequency combs emerged in micro-mechanical resonators [4,5] due to their strong nonlinearity [6,7].
Nonlinearity can be also present when distinct particles and quasiparticles [8] interact. This has led to the realization of frequency combs utilizing, e.g., acoustic and superconducting resonators [9] and magnon-phonon interactions with a Yttrium Iron Garnet (YIG) sphere in a cavity [10]. Magnetic materials are intrinsically nonlinear [11,12] so that the quanta of angular momentum, the magnons, have been shown to exhibit frequency combs when driven to nonlinear interactions [13], when interacting with solitons [14], driven by two tones [15], or via magnetostriction [16].
Another class of materials where frequency combs could be realized are articial spin ices [17]. These are geometrically arranged nanomagnets that host a wealth of dynamical modes [1820] including nonlinear dynamics [21].
Many other congurations are possible [17],
including extension to three-dimensional lattices [22],
which have also shown interesting dynamics such as coherent waves in tripod-based ASIs [23,24] and ultrastrong magnon coupling in 3D square ASIs [25]. How- * dbozhko@uccs.edu † eiacocca@uccs.edu ever, in all these cases, nanoscopic arrangements are spatially xed due to the lithography processes used, and the GHz dynamics are coupled through edge modes in which the magnetic volume and the concomitant dynamic stray eld is greatly reduced in extent. Therefore, frequency combs are likely unavailable in ASIs unless additional degrees of freedom are allowed.
Here, we demonstrate that a frequency comb can emerge from a driven macroscopic mechano-magnetic articial spin ice system [2628], or macro-ASI. By using a macro-ASI, we ensure a strong coupling mediated by stray elds throughout the array. In addition, the mechanical motion of the magnets results in a temporal modulation of the dipole-dipole interaction. When driven to nonlinear resonance, the system enters an unstable region characterized by a nonlinear amplitude and phase modulation that establishes the comb spectrum. A notable feature of our macro-ASI system is the richness of nonlinear dynamics at macroscopic scales and frequencies perceivable by the naked eye. Our results suggest that frequency combs can be also achieved in geometrically positioned arrays of micro-mechanical resonators with high quality factors (e.g. [4,29]) decorated with magnetic materials.
We used the square ASI conguration [30] shown in (b) (c) 10.4 Hz 10.4 Hz 11.8 Hz 11.8 Hz 8.68 Hz 8.46 Hz (a) d FIG. 1. (a) Picture of the macro-ASI composed of 60 magnets arranged in a square lattice. Resonant modes obtained from (b) experiments and (c) numerical modeling. The modes are in good quantitative agreement, both in their spatial symmetries and resonant frequencies. The scale bar represents 5 cm. ASI fabrication are given in the Appendix A. In contrast to previous realization of macroscopic ASIs [26, 27, 32], our set-up was built with a focus on driven dynamics rather than population of defects due to frustration from a high-energy state. Additional attention was paid to achieve an accurate lattice geometry and uniformity of the mechanical properties of individual elements. As a result, the macro-ASI exhibits wave propagation due to a localized perturbation, as shown in the Supplementary Video 1 [33]. Because of its nite size, the system supports only discrete standing wave modes of oscillation. These are experimentally measured by perturbing the ground state with a magnetic eld that approximates a Heaviside function produced by a solenoid of radius a = 10.5 cm placed directly under the macro-ASI. The time-dependent dynamics of each magnet's rotation angle φ was optically detected by a high-speed camera with a frame rate of 250 frames per second after the eld is rapidly turned o. The experimental details on the excitation as well as the data acquisition and processing is described in the Appendix B. Fourier analysis is performed on each time trace to obtain the frequency-dependent amplitude as a function of space. The reconstructed standing wave patterns are shown in Fig. 1(b), where the Fourier peak amplitude for each magnet was used to scale a twodimensional Gaussian function for better visualization. The symmetry of the system is apparent from the resulting standing wave patterns. The modes at 8.68 Hz and 10.4 Hz both exhibit excitation at the center of the array and symmetric excitation of magnets towards the diagonal edges of the array. The mode at 11.8 Hz corresponds to a bulk mode.
These standing wave patterns are well reproduced by a torque equation that considers the torque on each magnet due to the dipole-dipole interaction between all magnets as well as mechanical friction and inertia [26,27]. The equation reads
where I is the moment of inertia and φ i is the angle of monopole i relative to the z-axis. The rst term in the right-hand side describes the torque due to a the magnetic ux density B(x, y, z) that, in general, is a function of space. We dene L μi as the distance from the rotor to monopole i, which is free to rotate in the μi - k plane. The second term in the right-hand side describes the coupling between the monopoles i and j throughout the N magnets in the array. Interactions between the same monopoles in the same magnet are not allowed. The distance between monopoles is r μi,μj = L μj + d d -L μi , where d l is the distance between the rotors of magnet i and j. This implies that the distance varies in time as the magnets deviate from their equilibrium direction, providing a nonlinear coupling between the mechanical and magnetic forces. In general, the eective torque on the magnet is a three-dimensional vector. However, the magnets are constrained by the rotors so that only the μi × k component of the torque leads to motion. Note that this cross product also ensures an appropriate denition of the torque sign in the Cartesian reference frame for magnets oriented along either the x or y axes. Finally, we also include energy dissipation through friction scaled by the parameter η. We estimate η = 10 -7 Kg m 2 /s for a relaxation of approximately 2 seconds. Fur- ther discussion on the eect of the friction parameters is found in Appendix C. Equation (1) was solved using MATLAB's ode15s useful for sti dierential equations. We set outputs from the solver at increments of 0.01 seconds and ran it for 50 s to obtain a frequency resolution of 0.02 Hz and a maximum frequency of 50 Hz. Further details on our implementation of the model is given in the Appendix D. The resulting standing wave patterns from this model are shown in Fig. 1(c). To achieve this excellent agreement, we tuned the saturation magnetization of the magnets to M s = 899 kA/m, in good agreement with the 863 ± 15 kA/m found by SQUID magnetometry (detailed in the Appendix E). The remaining physical parameters were obtained from the nominal dimensions of the magnet: diameter D = 1.59 mm, length 2L = 25.4 mm, and mass 337 × 10 -6 kg. From these parameters, we calculate the magnet's moment of inertia for a cylinder I = 2.03 × 10 -8 kg m 2 . III. DRIVEN DYNAMICS AND FREQUENCY COMB The harmonically-driven dynamics of the macro-ASI were achieved by passing an alternating current produced by a function generator through the solenoid, as described in the Methods. The eld-magnitude-dependent spectrum at a drive frequency of f exc = 10.2 Hz is shown in Fig. 2(a), exemplary for one magnet of the array. At low ux density magnitudes, only the direct excitation is observed, indicating a linear dynamical regime. However, a change is observed at the threshold value of B 1 = 0.104 mT, above which the spectrum displays equidistant harmonics. This feature indicates that the dynamics become nonlinear, resulting in the emergence of a frequency comb. For example, the spectrum at B = 0.112 mT is shown in Fig. 2(b) for an extended frequency range. Up to thirteen peaks can be clearly discerned between 5 Hz and 12 Hz. Note, that the optical measurement method utilized in this study is linear with respect to the magnet's angle and does not distort the measured dynamics as an electromagnetic pickup coil would. The onset of the frequency comb depends also on the drive frequency. The frequency-dependent spectrum at our maximum ux density of B = 0.13 mT is shown in Fig. 2(c). The comb is only visible in the vicinity of f exc ≈ 10 Hz. Over this drive frequency, harmonic content about f exc /2 appears with well-dened features. The appearance of these half-frequency excitations in the spectrum could be understood as a result of a parametric instability of the forced oscillations into a pair of standing wave modes, which exist in our system at the half of the drive frequency. Further investigations of this interesting phenomenon requires more detailed insight into the dispersion relation of propagating waves in this system and is outside of the scope of this particular work. We note that these spectra have a striking resemblance to those obtained numerically for a nanoscopic square ASI driven to nonlinear regime [21], underscoring the relevance of mode mixing. However, in nanoscopic ASIs, the arrays are extensive and modes are located within the nanomagnets whereas in the present macro-ASI the modes are standing waves resulting from the constrained dimensions of the array. To gain further understanding of the behavior of the system, we performed numerical modeling, which quantitatively reproduced the experimentally observed driven dynamics. To capture possible eects emerging from non-uniformity of the excitation eld, we considered the spatially-inhomogeneous eld distribution due to a solenoid aligned with the geometric center of the macro-
We investigate the eect of structural parameters in the occurrence of the frequency comb. To restrict our numerical parameter space, we consider the solenoid geometry and dimension as well as the eld magnitude of for larger arrays at roughly the same frequency but with a dierent range. This can be understood by the fact that the solenoid remains unchanged and so more energy must be given to the system in order to excite the whole lattice into nonlinear dynamics. For the small array, the comb is not observed. This is due to the fact that there are insucient degrees of freedom to access the modes.
The bottom row of Fig. 4 shows the dependence on the lattice constant. In this case, the changes are more dramatic. For closely-spaced magnets, the comb is absent because of the strong coupling again arresting nonlinear dynamics at the chosen eld magnitude. As the lattice constant increases, the coupling between the magnets is reduced, resulting in a general redshift of the resonant modes. When the shift is too large, the comb is not well dened, e.g., at d = 6.31 cm.
It is also apparent that the frequency comb spacing is reduced by increasing the lattice constant. If we consider that the comb requires a low-frequency mode to exist, then the general redshift would also apply to such a mode and thus the reduction of the frequency comb spacing. This is more clearly seen in Fig. 5 where we compare the spectrum of the frequency comb at d = 5.05 cm driven by 10.2 Hz and d = 5.68 cm driven at 7.4 Hz. The frequency axis was centered to easily compare the combs. There is a clear redshift in the frequency comb spacing estimated to be a modest 6 mH. Based on these simulations, it may be possible to further reduce the frequency spacing of the comb in a larger array driven by a larger solenoid and likely a stronger eld. Amplitude modulation is not sucient to produce a frequency comb. Therefore, we additionally investigated the relative phase between the magnet and the drive frequency from the modeling. This is achieved by computing the phase dierence between the drive signal b and the angle φ [34] When the eld decreases in magnitude, shown by red circles, the comb is absent. We now understand the comb as a result of a Hopf bifurcation in the system: the xed point solution transitions into a periodic orbit. Physically, the bifurcation occurs because the drive frequency excites coexisting modes with dierent phases. As a result, the magnets periodically transition from one mode to the other. The modes involved in this process are observed in Fig. 7. The panels show the mode proles at the (a) rst lower harmonic, (b) fundamental mode, and (c) rst higher harmonic. The lower harmonic exhibits a mode prole consistent with the 10.4 Hz natural mode.
This is reasonable because it is the closest natural mode to the 10.2 Hz drive frequency. The rst higher harmonic, however, is a uniform mode. This uniform mode is excited by the drive frequency so it is not observed in Fig. 1(b) or (c). The fundamental mode is a combination of these two modes. It is also worth noting that the mode is not mirror-symmetric. This broken symmetry allows for the modes to be metastable. This is similar to the process observed in a nanomechanical resonator [4] in terms of nonlinear resonance and also as frequency pulling [35]. However, our system is driven into nonlinearity due to the strong dynamic coupling between the magnets. In fact, the model shows that at elds over B ≈ 1.02 mT, linearity is restored due to the dynamic mode coupling becoming negligible.
The mode metastabiliy can be observed experimentally by initiating the dynamics without any particular protocol, as evident from the Supplementary Videos 2 and 3 [33].
Our ndings are in stark qualitative contrast to the magnetization dynamics of nanomagnets, which are afxed to a surface. Strong static coupling is present in ASIs, as evidenced by the rich literature on frustration, e.g., Refs. [17,30,36,37], and dynamics due to emergent monopoles [38,39]. However, the dynamic coupling allowing for the wave propagation is extremely weak because it is due to edge modes which have a small mode volume [18,4042]. In fact, the edge mode volume cannot exceed a third of the nanomagnet's volume [42] and it is typically found in micromagnetic simulations to be much more localized [21]. As a consequence, the bulk modes are essentially decoupled as evidenced by their near-ideal Kittel dependence [43,44], even when forced into the nonlinear regime [21] or nanomagnets are stacked into 3D
structures [25]. This implies that the dynamic coupling in a nanoscopic ASI is approximately an order of magnitude weaker than that of a macro-ASI where the dynamic coupling is enabled by mechanical motion so that the static stray eld is modulated. Our results suggest the possibility to achieve richer dynamics in nanomagnets patterned on micro-mechanical oscillators, enabling dynamic coupling through the magnet's strong static stray eld and its modulation through high quality factor mechanical motion. Such devices would achieve a nonlinear regime at smaller driving elds due to the strong magnetic coupling between elements, possibly even in a quantum regime. From a broader perspective, our results demonstrate the generality of nonlinear phenomena giving rise to frequency combs, including macroscopic scales.
using precision Genmitsu 3018-PRO CNC desktop mill (improved in-house to achieve ≈ 10 µm accuracy).
Appendix B: Experimental data acquisition and processing
The experiment was set up by driving forced oscillations using an external magnetic eld. To achieve this, we introduced a magnetic coil with an outer diameter of The frequency of the driving harmonic voltage was systematically swept from 2.0 Hz to 20.0 Hz. To ensure the system was stabilized, we used to following protocol: rst, the frequency was linearly swept from 0 Hz to the desired measurement frequency in 30 s. Then the system was let to stabilize for 10 s. Finally, the measurements were taken by capturing a 30 s AVI video with the camera.
After the completion of the recording of the entire frequency sweep, the data was processed using a specially developed Python script based on OpenCV library [46].
Firstly, bounding rectangles containing each of the 30 magnets were cropped individually and then the blue and the red binary color masks were applied. Then, a line weighted t connecting the pixels was applied to estimate the magnet orientation. The deviation of the tted line from the horizontal axis indicated the angle of the magnet per frame. Each of these traces was analysed by fast Fourier transform (FFT). FIG. 10. The saturation magnetization of the magnets was tuned in the model to match the highest-frequency mode. The magnitudes of the spectral peaks are normalized to coincide at the 12 Hz mode. Excellent agreement is observed with modes at ≈ 8.5 Hz and 10.4 Hz. The experiment shows a much richer spectrum in the range 6 Hz to 10 Hz. The simulations overestimate the low frequency mode, due to the ideal placement and conditions at the boundaries.
the measured radius of a = 10.5 cm. We also consider that the solenoid is aligned with the geometrical center of the macro-ASI. The ux density amplitude B 0 is a controllable parameter in the modeling. We calibrate the eld by normalizing Eqs. (D1) at the centre of the solenoid and oset z. This means that r = z. Numerical computation results in a normalization factor of 0.8895.