Microwave resonances of magnetic skyrmions in thin film multilayers

Non-collinear magnets exhibit a rich array of dynamic properties at microwave frequencies. They can host nanometre-scale topological textures known as skyrmions, whose spin resonances are expected to be highly sensitive to their local magnetic environment. Here, we report a magnetic resonance study of an [Ir/Fe/Co/Pt] multilayer hosting Néel skyrmions at room temperature. Experiments reveal two distinct resonances of the skyrmion phase during in-plane ac excitation, with frequencies between 6–12 GHz. Complementary micromagnetic simulations indicate that the net magnetic dipole moment rotates counterclockwise (CCW) during both resonances. The magnon probability distribution for the lower-frequency resonance is localised within isolated skyrmions, unlike the higher-frequency mode which principally originates from areas between skyrmions. However, the properties of both modes depend sensitively on the out-of-plane dipolar coupling, which is controlled via the ferromagnetic layer spacing in our heterostructures. The gyrations of stable isolated skyrmions reported in this room temperature study encourage the development of new material platforms and applications based on skyrmion resonances. Moreover, our material architecture enables the resonance spectra to be tuned, thus extending the functionality of such applications over a broadband frequency range.

M agnetic skyrmions are stabilised due to competition between various spin interactions including Heisenberg exchange, the Dzyaloshinskii-Moriya interaction (DMI) and geometric frustration. These non-collinear configurations of magnetic moments have been envisioned as a platform for the next generation of spintronic devices. An essential requirement to realise this technological ambition is a complete understanding of the static and dynamic properties of skyrmions.
The majority of existing literature on skyrmion dynamics is restricted to non-centrosymmetric single crystals hosting either Bloch or Néel skyrmion lattices at low temperatures. These are magnetic materials with space group P2 1 3 (e.g. Cu 2 OSeO 3 , MnSi or Fe 1−x Co x Si), or R3m (e.g. GaV 4 S 8 ) 1,2 . In both types of system, the skyrmions exhibit clockwise (CW) and counterclockwise (CCW) gyration modes as well as a breathing (BR) mode [3][4][5][6][7][8] . In contrast with magnetic skyrmion arrays in bulk single crystals, Néel skyrmions at ferromagnet-heavy-metal (FM-HM) interfaces are stabilised by the interfacial Dzyaloshinskii-Moriya interaction. Their configuration can range from dilute to dense disordered arrays, as well as ordered lattices. However, the skyrmion size, stability and configurations are all sensitive to the magnetic field history and uniformity of the interface. Furthermore, additional energy contributions (such as long-range magnetostatic interaction between separate ferromagnetic layers) can also play a decisive role in determining the skyrmion properties, by modifying the helicity 9 . Magnetic multilayers are particularly promising candidates for developing skyrmion-based technologies, since they combine skyrmion stability at room temperature with the ability to finely tune their magnetic parameters via the multilayer geometry. The latter significantly amplifies the role of the dipole-dipole interaction, enriching the properties of skyrmion phases.
Recent experiments on FM-HM multilayers have yielded considerable insight into the static properties of Néel skyrmions [10][11][12][13][14] and their motion due to dc electrical currents [15][16][17] or magnetic field gradients 18 . However, the dynamic properties of Néel skyrmions in FM-HM heterostructures have so far resisted investigation, due to the reduced sensitivity of magnetic resonance probes in thin films and the difficulty of engineering high skyrmion densities. Consequently, many interesting and potentially useful predictions for skyrmion resonances 3,[19][20][21][22][23] , such as RF energy harvesters, microwave couplers and magnon gratings 2,24,25 , remain to be experimentally verified. A clear understanding of skyrmion responses at GHz frequencies is also prerequisite for operating recently-postulated synaptic computational architectures at useful clock frequencies [26][27][28] . Detailed experimental and numerical studies of skyrmion dynamics will be essential to deliver the many proposed applications requiring skyrmion interactions with radio or microwave frequency fields.
Here we present a broadband microwave absorption study of sputter-deposited multilayers hosting stable Néel skyrmions at room temperature and in out-of-plane (OP) fields 140 mT ≲ μ 0 H ⊥ ≲ 300 mT. Applying additional in-plane (IP) microwave fields in the range 5-15 GHz, we detect three resonant modes which correlate with the evolving magnetisation configurations in our samples. Our micromagnetic simulations indicate that the lowfrequency (LF) mode is associated with a CCW gyration of the skyrmion core, while the CCW precession of the magnetisation within the inter-skyrmion zones is the predominant contributor to the high-frequency (HF) mode. The intensity and resonant frequencies of both HF and LF modes are acutely sensitive to the interlayer dipolar interaction between ferromagnetic layers, which can be modulated by the thickness of the nonmagnetic spacer layer. The third mode emerges at high fields and is associated with the uniform precession of the field-polarised magnetisation.

Results
Our samples are [Ir 1 Fe 0.5 Co 0.5 Pt 1 ] 20 multilayers grown by magnetron sputtering, where the subscripts are the sequential layer thicknesses in nanometres, and the complete 'stack' comprises 20 repeats of [Ir 1 Fe 0.5 Co 0.5 Pt 1 ]. This stack is known to host Néel skyrmions from room temperature down to 5 K 14,29 . To probe the electrodynamic properties of these stacks, we employ a broadband technique [30][31][32] and measure the microwave transmission S 12 through a proximate coplanar waveguide (Fig. 1a) sequence of magnetic orders, revealed by magnetic force microscopy (MFM) imaging (Fig. 1b). For μ 0 H ⊥ ≳ 300 mT the magnetisation is saturated, but a slight decrease in field below~300 mT results in skyrmion nucleation at random sites. As we continue to reduce the field, further skyrmions are progressively nucleated, transforming the magnetic configuration into a densely-packed disordered skyrmion lattice which persists down to~140 mT. The disordered configuration most likely originates from local inhomogeneities in the magnetic interactions, due to structural defects and a complex energy landscape. Microwave resonances can be identified in all applied fields as local minima in the frequencydependent transmission S 12 (f) (Fig. 1c). Figure 2a depicts the microwave absorption across the phase diagram, measured via frequency-sweep spectroscopy at constant field. Searching for local minima in cuts along both field (Fig. 2b) and frequency (Fig. 2c) axes we identify three distinct resonant branches which correlate with the different magnetic configurations imaged via MFM. In the high-field polarised ferromagnetic state, the Kittel resonance displays the expected linear f(H ⊥ ) relation. Following skyrmion nucleation below~300 mT, the Kittel mode splits into two branches: a narrow resonance which softens rapidly, and a broad absorption which retains an approximately constant frequency as the field is reduced to~250 mT, below which the resonance frequency rises sharply with decreasing field. We label these branches LF and HF, respectively. For comparison, Fig. 2d shows absorption spectra from a similar heterostructure during field-sweep measurements. Despite the known sensitivity of skyrmion configurations to the magnetic field history 33 , the results are remarkably consistent, with all three resonances emerging at similar frequencies in both experiments. The LF mode cannot be resolved in linecuts along the frequency axis in either dataset (Fig. 2c, f), due to its relatively low absorption intensity and steep frequency dispersion with magnetic field. However, it is clearly visible in field axis linecuts (Fig. 2b, e) for data acquired in frequency and field-swept measurements. Our results therefore indicate the robust nature of both LF and HF resonances in the skyrmion phase.
To understand the microscopic origin of the experimentally observed LF and HF modes, we perform micromagnetic simulations using MuMax 3 34 . The magnetic configuration at each field was obtained by relaxing from a random magnetisation. Close to zero field, the magnetisation relaxes into a labyrinth stripe structure, which transforms into a dense skyrmion phase for μ 0 H ⊥ ≈ 140 mT. The evolution of the simulated OP magnetisation m z lies within 30 mT of our experimental results, as shown in Fig. 3a, confirming the accuracy of our simulation parameters. Slight deviations between the two datasets are to be expected since thermal fluctuations and spatial inhomogeneities (which may influence relaxation processes) are not considered in the present simulations. The spatial magnetic configuration is depicted in Fig. 3b. At intermediate fields even this clean system relaxes into a disordered skyrmion lattice, emphasising the complex energy landscape which can easily trap skyrmions in metastable configurations. With increasing field, the skyrmion array gradually transforms from a disordered lattice into a dilute gas of isolated skyrmions and, eventually, into a field-polarised magnetic state, in accordance with our MFM imaging (Fig. 1b). The validity of our simulation parameters-especially the exchange A and DMI D-is further confirmed by the close agreement of the simulated zero-field stripe periodicity~125 nm with the experimentally measured value~129 nm (Fig. 3c, d).
The simulated microwave absorption spectrum as a function of field and frequency is shown in Fig. 3e for an intrinsic damping parameter α = 0.05 determined from fits to our experimental data (Supplementary section I.A). Three distinct resonance modes can be identified in the simulated spectra, in agreement with experiments. The frequency sweeps indicate that the broad linewidth leads to a superposition of resonance signals, thus hampering a clear delineation of the LF and FM modes. For this reason, we performed further simulations with a lower damping parameter α = 0.01 (Fig. 3f), allowing us to resolve the dispersion of all the modes. We note that the simulated absorption intensity of the LF mode is higher than in our experimental data (Fig. 2). This can be attributed to the presence of extrinsic, modedependent damping mechanisms in our heterostructures which are not accounted for in our simulations (see Supplementary section I.B for further discussion). Nevertheless, the fielddependent resonant frequencies of all three modes are in close agreement with our experimental findings. As discussed below in more detail, the interlayer dipolar coupling is essential to reproduce our experimental results, emphasising the critical role of dipolar interactions in multilayers.
To elucidate the character of the skyrmionic excitation modes we performed further simulations on two additional configurations. In the first instance, a metastable state was prepared consisting of a single skyrmion in a field-polarised background. The absorption intensity of this specific state and its evolution with μ 0 H ⊥ is displayed in Fig. 4a. For the full range of applied fields 150 mT ≤ μ 0 H ⊥ ≤ 350 mT, the field-polarised background contributes to the spectral weight of the Kittel resonance, corresponding to the uniform precession of the field-polarised magnetic moments. At μ 0 H ⊥ = 280 mT, the size of the skyrmion suddenly increases with decreasing field, coinciding with the emergence of an additional resonant mode below the Kittel frequency. This resonance corresponds to a magnon-skyrmion bound state with CCW character (Fig. 4c) whose frequency tracks the LF resonance observed experimentally. Its spin-wave eigenfunctions extracted from the numerical data can be parametrised by the semi-major and semi-minor axes of the local precessional ellipse, f(r) and g(r), respectively, and agree well with analytical results (see Fig. 4b and Supplementary section III.C for details). Note that the sign of the product fg determines the sense of rotation of the magnetisation vector m around its equilibrium direction. Importantly, the spin wave functions exhibit a maximum close to the skyrmion radius but then decay rapidly to zero at large distances, indicating the localisation of the LF resonance within the skyrmion area. The corresponding time evolution is illustrated in Fig. 4c, where the large arrow represents the total dipolar moment performing CCW gyration.
Having identified the LF resonance as a CCW magnonskyrmion bound state, we now discuss the HF resonance. The frequency of the HF mode is located above the Kittel resonance frequency and thus within the scattering continuum of the fieldpolarised background, which suggests that it cannot be accounted for by the first setup containing only a single skyrmion. The HF resonance is therefore likely to be associated with a multiskyrmion configuration. We have verified that both HF and LF resonances observed in the simulations in Fig. 3 persist after annealing the metastable disordered skyrmion configuration into a regularly ordered hexagonal lattice. The annealing was achieved by shaking the system with the help of a large amplitude ac magnetic field, using a similar approach to vortex matter in superconductors 35 (see Supplementary section II.D).
In order to study the HF mode, we consider a second simulation configuration, consisting of a rectangular unit cell in a regular hexagonal lattice containing two skyrmions. Figure 5a shows the corresponding equilibrium magnetisation m 0 (x, y) at 200 mT. The numerically-obtained absorption spectrum of this system is depicted in Fig. 5b, together with the experimentally observed frequencies. In the field range encompassed by the skyrmion phase, two distinct resonances can be identified whose field dependences closely resemble those of the HF and LF modes observed experimentally. Moreover, the spectral intensity of the HF mode is stronger than that of the LF mode (Fig. 5c), in agreement with experiment. The magnon probability distribution of the two modes is obtained numerically by time-averaging hjmðx; y; tÞ À m o ðx; yÞj 2 i t and shown in Fig. 5d, g. The distribution of the LF mode is localised and concentrated on a ring around each skyrmion, in agreement with the result of the single-skyrmion calculations of Fig. 4b. In contrast, the probability density of the HF mode is mainly distributed within the area between the skyrmion positions. The precessional motion of the polarised magnetic moments in the inter-skyrmion zones also results in a CCW rotation of the total magnetic dipole moment for the HF mode (illustrated in our Supplementary videos).
A characteristic advantage of multilayers compared to noncentrosymmetric bulk crystals and thin-film monolayers is the tunability of the interlayer dipolar interaction by varying the stack geometry. To study the influence of interlayer dipolar interactions numerically, we varied the thickness L of the nonmagnetic spacer layer, corresponding to the [Ir-Pt] thickness in our multilayers. Thicker layers result in a relatively weaker interlayer dipolar interaction. The evolution of the resonance frequencies and absorption intensities as a function of L are shown in Fig. 5e, f, h, i. We find that the resonances are highly sensitive to the dipolar coupling when L is comparable to the thickness of the magnetic layer (1 nm), but saturate at a thickness L ≥ 10 nm, which is much smaller than the skyrmion diameter (~60-80 nm) or the inter-skyrmion distance (~130-150 nm). The resonant frequencies of the HF and LF modes respectively decrease and increase with L, before saturating at larger thickness. The spectral weight of the HF mode possesses a characteristic maximum at L~2-5 nm depending on the applied field, before decreasing and eventually vanishing for large L. The HF mode therefore cannot be observed in the limit of large L where the layers are practically decoupled. In contrast, the weight of the LF mode increases with L and saturates at a finite intensity for large thickness. Whereas the spectral intensity of both modes increases with applied field for thinner spacers, this trend reverses for larger L. For our experimental thickness L = 2 nm the intensity of the HF mode is larger than that of the LF mode for all fields, in agreement with our experimental results. This evident sensitivity of skyrmion resonances to the geometry of chiral magnetic multilayers highlights their potential for use in tunable microwave receivers.

Discussion
It is instructive to compare our results with previous work on bulk magnets with a Dzyaloshinskii-Moriya interaction [3][4][5][6][7][8] . The magnetic skyrmion lattice of these systems at low temperatures is known to possess a CCW mode and a clockwise (CW) mode for IP ac magnetic fields, as well as a breathing mode for OP ac fields. Such behaviour is distinct from the spectra we measure for skyrmions in multilayers: in particular, there is no CW mode within our experimental frequency range. This difference can be attributed to the presence of an uniaxial anisotropy K, which is absent in cubic chiral magnets, and-more importantly-to the strong influence of the dipolar interaction. The relative strength of the latter is quantified by the dimensionless parameter μ 0 AM 2 s =D 2 , where μ 0 is the vacuum permeability, A the exchange stiffness, D the DMI and M s the saturation magnetisation. In our multilayers this parameter is one order of magnitude larger than in bulk The colour scale and arrows represent the out-of-plane and in-plane magnetisations. b Simulated resonance spectra for the geometry described in (a). c Vertical linecut through the spectra in (b) at 200 mT, highlighting the relative intensity of the HF and LF modes. d, g Spatial distribution of the microwave absorption intensity during HF and LF resonances at 200 mT, excited at 10.6 and 2 GHz, respectively. While the distribution of the LF mode is localised close to the skyrmion core, the HF mode intensity is primarily concentrated in the inter-skyrmion area. Panels a-d and g are simulated with a spacer thickness L = 2 nm. Panels e, f and h, i depict the variation of the absorption intensities and resonant frequencies with L for the HF and LF modes, respectively. The red dashed line at L = 2 nm corresponds to the experimental layer separation. chiral magnets hosting Bloch skyrmions 2 . Moreover, Néel skyrmions possess finite volume magnetostatic charges, unlike Bloch skyrmions. These aspects drastically modify the excitation spectrum of skyrmion resonances for our magnetic multilayers in comparison with bulk chiral magnets. As a result, the CW mode has a significantly lower intensity and is shifted to high frequencies beyond our experimentally accessible frequency range. At the same time, the CCW excitation softens and is shifted to lower frequency due to this enhanced attraction of spin waves in the CCW channel. This produces the magnon-skyrmion bound state for a single skyrmion shown in Fig. 4, which is responsible for the LF resonance observed in our multilayers. A bound CCW state of this nature does not exist in cubic chiral magnets 21 . In addition, due to the enhanced attraction in the CCW channel, the next higher-order CCW mode (possessing an additional node in its wavefunction, as shown in Fig. 5d) becomes experimentally accessible, resulting in our observed HF mode. This pronounced sensitivity to dipolar interactions enables a remarkable tunability of the CCW resonances through the multilayer geometry. We note that our numerical calculations also predict a breathing mode in the multilayer with lower frequency~1.5 GHz. However, this can only be probed with OP ac magnetic fields (Supplementary section II.F) and is beyond the experimental scope of this study.
In summary, we have studied the dynamic properties of magnetic skyrmions in technologically relevant multilayers hosting Néel skyrmions at room temperature. Broadband microwave absorption measurements reveal the presence of a low-frequency resonance associated with CCW gyration of isolated magnetic skyrmions. The interlayer dipolar interaction in our multilayers generates a high-frequency collective skyrmion resonance with CCW character. Crucially, these resonant modes can be tuned over a broadband frequency range by minor adjustments in the stack geometry. Combined with our earlier results 10 , we have now established that both static and dynamic properties of skyrmions can be tuned via the multilayer stack architecture. Our results expose the potential of skyrmion-hosting multilayers as a material platform capable of combining conventional dc spintronics (for data storage and logic operations) with magnonic or synaptic devices operating at GHz frequencies.

Methods
Sample growth and initial characterisation. All our [Ir/Fe/Co/Pt] heterostructures are grown on Si wafers, with a~1 μm-thick thermally oxidised surface layer. The multilayers are grown by dc magnetron sputtering in an ultra-high vacuum chamber. The multilayer composition is Ta 3 /Pt 10 /[Ir 1 /Fe 0.5 /Co 0.5 /Pt 1 ] 20 / Pt 2 , where the subscripts indicate the layer thicknesses in nm and the magnetically active [Ir 1 /Fe 0.5 /Co 0.5 /Pt 1 ] unit is stacked 20 times. The Ta 3 /Pt 10 seed layer is included to optimise the crystallinity of the upper layers. Deposition at low power (25 W) and chamber pressure (1.5 × 10 −3 Torr Ar) allows us to control the layer thicknesses with Angstrom sensitivity, resulting in systematic and reproducible magnetic properties of our heterostructures. We determine the magnetic parameters of our multilayers using a combination of bulk magnetometry and local imaging. A Quantum Design TM MPMS-XL SQUID magnetometer was used to measure the saturation magnetisation M s . Spatially resolved magnetic imaging was performed using a Bruker TM D3100 atomic force microscope equipped with Nanosensors TM ultra-low moment magnetic tips (radius < 15 nm) and a homogeneous, adjustable perpendicular field provided by a permanent magnet array. Images were typically acquired with a tip lift height of 20 nm, which is sufficient to prevent the stray field of the tip from influencing the chiral spin textures.
Ferromagnetic resonance measurements. We acquire field-sweep and frequency-sweep microwave absorption spectra via a broadband technique, previously used to measure resonance in a range of thin-film magnets. Our heterostructures are secured upside-down on a coplanar waveguide using a spring-loaded clamp. This waveguide is mounted on the cold finger of a variable-temperature probe in an electromagnet, and the microwave transmission S 12 is measured using a Keysight TM PNA N5222 vector network analyser (VNA) in two-port mode. The useful bandwidth of this technique (primarily limited by the length of the coaxial cables) extends above 30 GHz, far beyond the frequency range studied here. Fieldswept measurements exhibit lower noise (and a larger normalised drop in S 12 at resonance), due to the limited response time and hence reduced sensitivity of the VNA during frequency sweeps. The measured FMR spectra contain superposed absorptive and dispersive components, the latter contributing to the phase lag between the microwave excitation field and the magnetisation response 36 . Following Dyson's original approach to modelling spin resonances 37 and standard practices in the FMR community [38][39][40] , we fit all our experimental spectra with a superposition of symmetric and antisymmetric functions: where A is the amplitude of the (Lorentzian) absorption component, β is the dispersion to absorption ratio, D accounts for a linear drift in the VNA output over time 40 and C describes the constant background signal from cable losses. Resonance occurs at x = x 0 with spectral linewidth Δx, where x ≡ f for frequency sweeps and x ≡ H for field-swept data.
Micromagnetic simulations. We model the local magnetisation and microwave response of our multilayers using the MuMax 3 software package 34 . This software solves the Landau-Lifshitz equation ∂ t m ¼ À γ 1þα 2 ðm B eff þ αm ðm B eff ÞÞ for the unit magnetisation vector m = M/M s with gyromagnetic ratio γ, damping constant α and saturation magnetisation M s . The random stochastic field which can be used to simulate thermal fluctuations was switched off in order to avoid unnecessary noise. The effective magnetic field B eff = −δE/δM is derived from the magnetic energy E, which is explicitly shown in the supplementary materials (section III.A) and includes the exchange interaction A, DMI D, uniaxial anisotropy K, Zeeman energy and dipolar interaction. The parameters A = 9.25 pJ/m, D = 1.40 mJ/m 2 , K = 0.65 MJ/m 3 and M s = 1.02 MA/m were adjusted to provide the best fit to our combined experimental magnetometry, MFM and FMR data (see Supplementary section I.A). The discretisation is chosen in the form of a cuboid with dimensions Δx × Δy × Δz, where Δz is set to 1 nm for all geometries (except for the thickness dependence study in Fig. 5 where Δz = 0.5 nm for L = 0.5 nm). For the large 2 × 2 μm 2 geometry simulations of Fig. 3, the in-plane discretisation is Δx = Δy = 3.9 nm, which is reduced to 1 nm for the additional simulations in Figs. 4 and 5. We simulated a ferromagnetic-nonmagnetic [FM (1 nm)-NM (L)] bilayer with L = 2 nm corresponding to the experimental setup. The thickness L was also varied for the simulations presented in Fig. 5e, f, h, i. Periodic boundary conditions (PBC) were applied in all three directions, thus approximating the 20 bilayers of the experimental geometry. For every magnetic field step in our simulations, the system was relaxed from an initial random magnetisation state. This yields disordered skyrmion lattice configurations which are metastable, since they can be annealed into an ordered skyrmion lattice by the repeated application of oscillating magnetic field pulses. As these metastable configurations resemble our experimental observations, they have been used to obtain the net magnetic moment, the magnetisation maps and the microwave absorption spectra shown in Fig. 3. For the spectra, we applied a sinc pulse of an excitation magnetic field b(t) = e x b 0 sinc(2πf max (t − t 0 )) with f max = 20 GHz, t 0 = 1 ns and b 0 = 10 mT, then Fourier transformed the induced in-plane magnetisation oscillations m x . For the sake of comparability with the experimental data, we show the absolute value of the Fourier transform. Noise at small frequencies is tentatively attributed to the relaxation processes of the metastable state. The simulations in Fig. 4 were performed on an in-plane 1 × 1 μm 2 area with a metastable configuration comprising a single skyrmion in a field-polarised background; a sinc pulse with b 0 = 2 mT was applied to generate the absorption spectra. Finally, the simulations of Fig. 5 were performed on an in-plane area corresponding to a unit cell of a hexagonal skyrmion lattice, where the lattice constant was adjusted as described in the main text. Here, a sinc pulse with b 0 = 1 mT was applied to obtain the absorption spectra.

Data availability
Data are available from the corresponding authors upon reasonable request.