Scanning Spin Probe Based on Magnonic Vortex Quantum Cavities

Performing nanoscale scanning electron paramagnetic resonance (EPR) requires three essential ingredients: First, a static magnetic field together with field gradients to Zeeman split the electronic energy levels with spatial resolution; second, a radio frequency (rf) magnetic field capable of inducing spin transitions; finally, a sensitive detection method to quantify the energy absorbed by spins. This is usually achieved by combining externally applied magnetic fields with inductive coils or cavities, fluorescent defects, or scanning probes. Here, we theoretically propose the realization of an EPR scanning sensor merging all three characteristics into a single device: the vortex core stabilized in ferromagnetic thin-film discs. On one hand, the vortex ground state generates a significant static magnetic field and field gradients. On the other hand, the precessional motion of the vortex core around its equilibrium position produces a circularly polarized oscillating magnetic field, which is enough to produce spin transitions. Finally, the spin–magnon coupling broadens the vortex gyrotropic frequency, suggesting a direct measure of the presence of unpaired electrons. Moreover, the vortex core can be displaced by simply using external magnetic fields of a few mT, enabling EPR scanning microscopy with large spatial resolution. Our numerical simulations show that, by using low damping magnets, it is theoretically possible to detect single spins located on the disc’s surface. Vortex nanocavities could also attain strong coupling to individual spin molecular qubits with potential applications to mediate qubit–qubit interactions or to implement qubit readout protocols.


I. INTRODUCTION
Electron Paramagnetic Resonance (EPR) is widely used in chemistry, physics, medicine and material science to characterize the electronic structure of magnetic molecules and impurities 1 .This has important applications in the study of organic and inorganic free radicals, colored centers in crystals, tissue oxygenation and archaeological dating, to give a few examples.Similarly to the physics of nuclear magnetic resonance imaging 2 , EPR can be combined with nanoscopic field gradients to detect spatially distributed spins 3 .
The technological interest in EPR has led to the development of very sophisticated and sensitive detection methods.These include optical techniques, such as using optically active atomic defects [4][5][6][7] , electrical detection of magnetic signals using scanning tunneling microscopy probes [8][9][10] or mechanical sensing based on magnetic resonance force microscopy 11,12 .Among the inductive methods, pickup coils and superconducting quantum interference devices (SQUIDs) have been employed to characterize resonant phenomena in small paramagnetic crystals 13,14 .However, the most widespread inductive-EPR readout technique uses transmission lines and cavities [15][16][17] .This latter approach offers the advantage of confining light in the space domain, yielding ina) Electronic mail: dzueco@unizar.esb) Electronic mail: pemar@unizar.escreased light-matter interactions.In addition, the cavity's quality factor can be further enhanced when using superconducting coplanar resonators instead of metallic three dimensional cavities 18 , yielding increased visibility.This idea is exploited in circuit Quantum Electrodynamics (QED) 19 , allowing, e.g., the manipulation and interrogation of superconducting or magnetic qubits 20 and quantum sensing of small amounts of spins [21][22][23][24][25][26][27][28][29] .
The light-matter coupling factor (g) depends on the amplitude of the (position-dependent) root mean square (rms) vacuum magnetic field fluctuations in the cavity B rms (r).The latter increases for decreasing mode volume, favouring the use of small cavities.However, downsizing is limited by the maximum operating frequency (1 − 20 GHz, typically) and the impedance of the resonator.Large coupling factors can be reached by fabricating nanoscopic constrictions at the central transmission line in superconducting coplanar resonators [30][31][32] .This approach keeps the total length of the cavity while reducing the other two dimensions, yielding strong focusing of the rms vacuum field fluctuations in nanometric regions around the central conductor.By doing so, large coupling factors of g/2π ∼ 1 kHz per spin have been demonstrated 33 .Alternatively, low impedance LCresonators allow increasing the amplitude of the magnetic compared to the electric field fluctuations, also yielding increased couplings 22,23,26,29,34 .
In addition to photons, the solid state offers a wide variety of bosonic excitations that can be emitted or absorbed such as, e.g., quantized spin waves or magnons [35][36][37][38] .In what follows, we will use g to de-note both the spin-photon and the spin-magnon coupling.Magnonic cavities could be used to perform spin qubit readout or to mediate spin-spin interactions [39][40][41][42][43][44][45][46][47][48][49] , offering the advantage of increasing the coupling by operating at reduced wavelengths (compared to electromagnetic resonators of the same frequency).This is possible since spin wave modulation is only limited by the lattice constant of the ferromagnet, allowing the downsizing to the nanometer range.This significantly decreases the mode volume and results in substantial spin-magnon couplings.For example, quasi-homogeneous spin waves in saturated ferromagnets have been proposed to substitute superconducting cavities 50,51 .Using nanoscopic Yttrium-Iron-Garnet (YIG) spheres, such approach shall provide strong coupling to individual free electrons, even reaching g/2π ∼ 1 MHz.Interestingly, whispering gallery modes in relatively large vanadium tetracyanoethylene (V[TCNE] x ) discs would yield a sizable spin-magnon coupling to individual NV centers 52 .This is so for the relevant mode volume is given by the disc's thickness t and the angular index magnon mode.In this way, a very encouraging g/2π ∼ 10 kHz can be obtained with a R = 500 nm, t = 100 nm out-of-plane saturated disc at 1.3 GHz.However, both approaches do require an externally applied bias field B ap with a double purpose: to saturate the ferromagnetic volume and to tune its resonance frequency with that of the spins.Interestingly, magnons can be confined in peculiar flux-closure configurations at B ap = 0.This is the case of magnetic vortices that are easily stabilized in thin-film ferromagnetic structures with lateral size between a few 10 nm up to several micrometers 53,54 .Minimization of magnetostatic energy yields a circular in-plane arrangement of spins with a nanoscopic out-of-plane magnetization core in the center.The vortex core modifies the spin-wave spectrum of the ferromagnet, adding new resonant modes in the absence of externally applied fields 55,56 .
Here, we compare the spin coupling to microwave resonators with that resulting from magnonic cavities (both saturated ferromagnets and flux-closure states).In the case of magnons, we base our calculations on three archetypical materials, common in the field of quantum magnonics [57][58][59] : On the one hand, a ferrimagnetic garnet having record low damping but low saturation magnetization (YIG, with Gilbert damping parameter α ∼ 10 −4 ).On the other hand, two ferromagnets with larger magnetization but higher damping, i.e., Co 25 Fe 75 alloy (CoFe, with α ∼ 10 −3 ) and Permalloy (Py, with α ∼ 0.5 × 10 −2 ), the most widely used soft-ferromagnet for spin-wave-based applications.We first demonstrate that the coupling to magnons is more than two orders of magnitude larger than that resulting from cavity photons.Although being of the same order, using vortices has important advantages over saturated ferromagnets such as the absence of a biasing magnetic field, independence from the saturation magnetization (M sat ), and the capability to manipulate the vortex's position.This ability enables scanning over a range of tens of nanometers.
In summary, our findings highlight the potential of fluxclosure states compared to other systems, demonstrating that the vortex core can be operated as a nanoscopic scanning EPR probe.Coupling the ferromagnetic disc to a superconducting circuit, our approach shall allow to spatially resolve the location of single spins distributed over the surface of the disc.

II. RESULTS
EPR is a resonant phenomenon, and regardless of the use of superconducting or ferromagnetic cavities, it must satisfy two conditions.The first criterion is energetic: the frequency ω 0 of the rf field produced by the cavity needs to match the energy difference between Zeemansplit spin levels.In the case of free S = 1/2 spins, the latter means ω s = γ e B tot (r) = ω 0 .Here, γ e /2π = 28 GHz/T and B tot (r) ≡ |B tot (r)| is the total (positiondependent) magnetic field that contributes to the Zeeman splitting.This criterion leads us to the definition of the resonance window, i.e., the region in space where spins are resonant with the cavity (see Section Methods).We highlight that the resonance window will only depend on the modulus of the (position-dependent) B tot (r).
The second condition is imposed by the geometry of the experiment: spin transitions can be induced only by the components of the vacuum field fluctuations that are perpendicular to the quantization axis of the spin.In the case of free S = 1/2 spins, the latter is parallel to B tot .This determines the strength of the spin-photon coupling g.To quantify the coupling, it is convenient to introduce the set of mutually orthogonal vectors B rms,j , j = 1, 2, 3 so that B rms,3 ∥ B tot .Thus, the other two components can induce spin transitions, allowing us to define B rms = B 2 rms,1 + B 2 rms,2 .Consequently, in the case of S = 1/2, the spin-photon coupling is given by g = µ B B rms with µ B being the Bohr magneton. 60We highlight that g will depend on both the intensity of the vacuum field fluctuations and also on the positiondependent distribution of its components with respect to B tot (r).

A. Coupling to spins, comparison between photons and magnons
We start by analyzing the spin-photon coupling in electromagnetic cavities.This will help us to compare with the case of magnons, which will be discussed below.We focus on the particular case of a superconducting coplanar waveguide formed by interrupting an open transmission line by two gap capacitors as sketched in Fig. 1a.The external magnetic field B ap = B tot shall be ideally applied along the transmission line axis so that, in this case, B rms = |B rms |.Under these circumstances, all spins can satisfy the resonance condition and are suscep- tible to be detected.In this case, the volume V of the resonator is the total length l (see Fig. 1a) multiplied by the cross-section of the central transmission line.l is fixed by the operating frequency ω 0 which, together with the circuit's impedance Z 0 , sets the intensity of the zeropoint current fluctuations 61 i rms = ω 0 ℏπ/4Z 0 .Using electromagnetic waves imposes a lower limit on the total length of the cavity l > λ/2 = 2πv/ω 0 , where λ and v are the wavelength and phase velocity, respectively.Working at high frequencies allows increasing the coupling by reducing the volume, but this is limited to ω 0 /2π < 10 − 15 GHz due to technical reasons.Decreasing V can only be accomplished by decreasing the cross-section which has the effect of focusing the vacuum field fluctuations while keeping i rms unchanged.Fig. 1b shows the spatial distribution of B rms for a ω 0 /2π = 10 GHz resonator of thickness 50 nm and different widths (5 µm, 500 nm and 50 nm).Fig. 1c shows B rms vs. V calculated at 3 nm above the central conductor for different values of ω 0 and the cross-section (10 × 10 nm 2 , 20 × 20 nm 2 , 35 × 35 nm 2 and 70 × 70 nm 2 ).From our simulations we see that the spin-photon coupling between free spins and Z 0 = 50 Ωresonators is limited to the range 1−50 kHz, in agreement with recent experiments 33 .Alternatively, spin waves propagate at much lower velocities, allowing to decrease λ while working at frequencies in the 1 -20 GHz range.This is the basic idea behind the use of magnonic cavities.One basic example to start with is the quasi-homogeneous Kittel mode (k = 2π/λ = 0) excited in isotropic ferromagnets, e.g., spheres.Zero point magnetization fluctuations in the ferromagnet produce zero point field fluctuations in the outside.Fig. 1e shows the spatial distribution of B rms for a saturated Py sphere for which impressive values of 200 µT can be reached.Fig. 1f shows B rms at 3 nm above the surface of spheres of different sizes (R = 25, 50 and 100 nm) and materials (CoFe, Py and YIG).The intensity of field fluctuations does not depend on B ap and satisfies sat , with V the volume of the ferromagnet.As it can be seen, the resulting couplings g = µ B B rms will be more than two orders of magnitude larger than those achievable with superconducting resonators, even approaching the 10 MHz range.
As a drawback, this approach requires unavoidably the use of an external bias field B ap that serves to saturate the ferromagnet and to tune its resonance frequency ω K = γ e B ap .In addition, unlike the case of a superconducting resonator, paramagnetic spins located close to the sphere's surface will feel a strongly non-homogeneous magnetic field B tot = B stray +B ap with B stray the demagnetizing field created by the sphere itself (see arrows in Fig. 1d).The latter means that not all spins will satisfy the resonance condition ω s = γ e B tot = ω K .As a matter of fact, spins located at the position of the white dot in Fig. 1e (where the coupling is maximum) will be far from resonance, therefore not contributing to the total signal.Independently of the magnitude of B ap , the resonance window for free paramagnetic spins will be fixed and it corresponds to the yellow volumes shown in Fig. 1d.
Much more interesting for applications is the study of flux-closure topological magnetic textures like the magnetic vortex sketched in Fig. 1g.Without requiring the application of external magnetic fields, vortex dynamics include the gyrotropic precession of the vortex core around its equilibrium position at frequencies in the range ω v0 ∼ 0.1 − 2 GHz.As an example, Fig. 1h shows B rms created by a R = 50 nm, t = 20 nm Py disc.We obtain very large B rms close to 200 µT as in the case of the Py sphere of similar volume (see Fig. 1e).Finally, we compute B rms at 3 nm from the disc surface.Results obtained for discs of different sizes (R = 50, 100 and 400 nm and thicknesses given in the legend) and materials (CoFe, Py and YIG) are shown in Fig. 1i.The intensity of the zero point field fluctuations is comparable to that obtained with the saturated spheres, also following the expected V −1/2 dependence.Unlike the previous case, B rms does not noticeably depend on M sat .In the following, we will analyze three distinct features that make vortex modes more suitable for spin detection compared to homogeneous modes.

i) Absence of biasing field: the resonance window
Vortex modes do resonate in the absence of any external biasing field B ap = 0, which is a great advantage compared to saturated ferromagnets.Besides, the demagnetizing stray field created by the vortex core itself reaches relatively large values close to the disc's surface (see dark arrows in Fig. 1g).This stray field is enough to make the vortex gyrotropic mode resonant with paramagnetic S = 1/2 spins fitting within the resonance window.In the case of the vortex, the latter takes the shape of a hollow semi-sphere (highlighted in yellow in Fig. 1g).The total effective coupling G will be enhanced by a factor √ N where N are the number of spins within the res-onance window (see Section Implementation of the vortex sensor).
The size of the resonance window depends on the radius of the vortex core itself r v .The latter is given by the material-dependent exchange length r v ≃ λ ex = 2A/µ 0 M 2 sat with A the exchange stiffness and µ 0 the vacuum magnetic permeability.r v is nearly constant, independent of the disc radius or thickness.However, for t ≫ λ ex , r v increases linearly with the thickness (see Fig. 2a).As a result, the resonance window increases considerably for materials with low M sat and large thickness.This can be seen in Fig. 2b, where we plot the spatial dependence of the spin resonance frequency ω s = γ e B tot for a paramagentic S = 1/2 spin as a function of its position above a R = 200 nm, t = 200 nm YIG disc and a R = 200 nm, t = 60 nm CoFe disc.The resonance window corresponding to the fundamental gyrotropic mode of the two discs is highlighted in red, i.e., the region for which ω s = ω v0 .As it can be seen, the resonance window of YIG is considerably larger (diameter ∼ 190 nm and height ∼ 28 nm) than that of CoFe (∼ 54 nm ×19 nm).This stems from the (more than one order of magnitude) smaller M sat of YIG compared to CoFe and the larger thickness of the YIG disc.

ii) Independence on Msat and higher order modes
The vortex core gyro, key to produce the B rms , also yields flexure oscillations of the vortex core line across the disc thickness 62,63 .For sufficiently thick discs, the latter leads to the emergence of higher order modes.This is shown in Fig. 2c where we plot the frequency spectrum of the same discs shown in panel (b).Apart from the gyrotropic fundamental mode (denoted v0), we find a first flexure mode (v1).At higher energies it is also easy to find azimuthal modes (a) that correspond to the rotation of two halves of the disc with opposite out-of-plane net magnetization.These modes can be distinguished in the the FFT of the spatially-resolved time-varying outof-plane magnetization in the upper and bottom surfaces and the transverse cut of the discs, as shown in the inset of Fig. 2c.Colored regions correspond to time-varying out-of-plane magnetization whereas darker areas correspond to constant magnetization.The signatures of vortex core precession can be seen in all modes.
Materials with low saturation magnetization are very interesting for sensing applications as lowest dampings are usually obtained in insulating magnets like YIG or V [TCNE] x .In this regard, the vortex mode offers an additional advantage over homogeneous Kittel modes for which g ∝ √ M sat .As shown in Fig. 1i, spin coupling to vortex modes is nearly independent of M sat encouraging the use of ultra low-damping ferrimagnets.However, low M sat yields vortex gyration at sub-GHz frequencies, usually too low for practical readout circuits.In this way, flexure or azimuthal modes shall allow operating ultralow damping ferrimagnets at reasonable frequencies.As an example, see the upper panel in Fig. 2b where we show the resonance windows (orange and yellow) for high frequency modes in a R = 200 nm, t = 200 nm YIG disc (ω v1 /2π = 760 MHz and ω a /2π = 1.7 GHz, respectively).Interestingly, the resulting spin-magnon coupling to high frequency modes is comparable to that of the fundamental mode.To illustrate this, we calculate the coupling for a single S = 1/2 spin located at the center, at 3 nm above the surface of the YIG disc.This yields g v0 /2π = 100 kHz, g v1 /2π = 99 kHz and g va /2π = 84 kHz for the fundamental, first flexure and first azimuthal modes, respectively.Additionally, high frequency modes in high M sat materials would resonate with paramagnetic spins lying much closer to the disc surface, where the coupling is larger.This reduces the total volume of the resonance window, boosting the spatial resolution of the vortex sensing probe.For example, in the case of the R = 200 nm, t = 60 nm CoFe disc, the resonance window reduces from ∼ 54 nm ×19 nm (mode v0) down to ∼ 31 nm ×9 nm (mode v1).

iii) Vortex mobility: Scanning spin detector
One of the most interesting properties of a spin sensor is the ability to scan over the surface of the sample.This is not trivial in the case of saturated ferromagnets but can be easily done in the case of magnetic vortices by using inplane magnetic fields.The effect of an in-plane external field B app∥ is to enlarge the disc region having magnetization pointing in the same direction while decreasing the size of the region having opposite magnetization (see inset in Fig. 3a).This results in an effective translation of the vortex core position d perpendicular to the direction of B app∥ .Fig. 3a shows the numerically calculated values of d vs. B app∥ for 60 nm-thick discs of different materials and radius.The vortex core moves progressively as B app∥ increases until it approaches the edge of the disc where it is annihilated.We highlight that the vortex state remains stable up to the annihilation field.Experimental measurements indeed demonstrate that the energy barrier for vortex annihilation remain substantial, reaching several hundred K, even at large magnetic fields of several tens of mT 54 .The latter holds true for discs of large radius but also for very small size, even below R = 100 nm.Additionally, the displacement of the vortex has little effect on the gyrotropic frequency as shown in Fig. 3b (left) 64 .The coupling strength and the size of the resonance window does not vary notably as long as the vortex displacement d < R/2.For larger displacements, the resonance window is considerably enlarged so that vortex fluctuations still allow spin detection.The resonance can be also modified by means of an out-ofplane magnetic field B app⊥ as shown in Fig. 3b (right) 65 .Interestingly, B app⊥ has also the effect of increasing (decreasing) the total size of the vortex core when applied parallel (anti-parallel) to the vortex polarity.Combining these two effects, it is possible to tune the total size of the resonance window.

B. Implementation of the vortex sensor
To finish, we discuss the experimental use of magnetic vortices to perform nanoscopic scanning imaging of a spin ensemble.The total spin-magnon coupling can be calculated as G = j g 2 j , where j goes from 1 to N spins satisfying the resonance condition.The sensitivity relies on the ratio between G and both the losses on the vortex resonator (γ v ) and the linewidth of paramagnetic spins (γ s ).Increasing the coupling is achieved by using discs of small volume since g ∝ V −1/2 (cf.Fig. 1i).However, the smaller the disc the more difficult its fabrication, manipulation and readout will be, so we will keep R ∼ 100−200 nm.The average coupling per spin defined as ⟨g⟩ = j g 2 j /N does not depend noticeably on the disc thickness (see Fig. 4a).For this reason, we will keep relatively large t ∼ 60 nm, favoring the stabilization of vortices with not too small resonance frequencies.
If the total effective spin-magnon coupling satisfies G > γ v , γ s , we are in the strong coupling regime and the resonance peak of the vortex mode splits into two peaks separated by 2G.Measuring the splitting provides, therefore, a direct way of detecting the presence of spins as, e.g., the vortex core is scanned.On the other hand, sensing will be also possible in the weak coupling regime.In this case, the width of the vortex resonance γ v will be enlarged if enough paramagnetic spins in the surface of the nanodisc satisfy the resonance condition.The linewidth γ v can be experimentally measured using superconducting microcircuits that can be optimally coupled to magnonic resonators 57,[66][67][68][69][70][71][72][73] .In this way, the target spins can be deposited in powder or crystal form or from solution on the surface of the disc for nanoscopic EPR imaging.
A suitable experimental approach is depicted in Fig. 4b.We consider a R = 200 nm, t = 60 nm Py disc that resonates at ω v0 /2π = 0.93 GHz with losses γ v = 14 MHz.We chose Py although being a lossy ferromagnet as it is probably the most common and easy material to fabricate and manipulate.The disc is located on the inductive part of a superconducting LC resonator of width 100 nm and thickness 50 nm, in good contact to it.The LC resonator is interrogated through a superconducting transmission line.Under such circumstances, the absorption dip produced by the resonant disc shall produce changes in the transmission that amount to several tens dB (see Fig. 4c).Now, we take the case of a typical EPR calibrating sample of free-radical molecules of 2-diphenyl-1-picrylhydrazyl (DPPH) having a density of ρ ∼ 2 spins/nm 2 , S = 1/2, g e = 2 and γ s /2π ∼ 10 MHz.Scanning the vortex sensor, one could detect the presence of small drops containing 2 × 10 3 DPPH spins (or, put in other words, volumes of only 0.4 atto-liter with 2 spins per nm 2 ), yielding G/2π = 1.2 MHz.The signature of coupling would be a variation of 0.1 dB of the resonance absorption as an external in-plane magnetic field is scanned (see Fig. 4b and inset in panel c).
Even more interesting for sensing applications, the use of low-damping materials like YIG opens the way to the detection of single spins.For instance, the first flexure mode in a R = 100 nm t = 60 nm YIG disc would bring one single spin (lying at 3 nm over the surface) into resonance at ω v1 /2π = 2.3 GHz, yielding G/2π = 0.7 MHz.For this purpose, an out-of-plane biasing field of 5 mT is necessary to tune the size of the resonance window.Un-der identical conditions as shown in Fig. 4b and thanks to the low losses of YIG (γ v = 0.6 MHz), this configuration would yield a large transmission signal depending on the relaxation time of the spin.This is exemplified in Fig. 4d where the transmission coefficient S 21 resulting for a bare disc and a disc coupled to one individual spin with γ s /2π = 10 MHz (0.3 dB difference) and 1 MHz (3 dB difference) are compared.

III. CONCLUSIONS
Our simulations suggest that, using Py, it is possible to detect small drops of 0.3 atto-liter containing 2 spins per nm 2 over a surface of 2π × 200 2 nm 2 .Additionally, the vortex core can be easily scanned by means of external magnetic fields so to perform EPR scanning microscopy.Interestingly, the same can be in principle achieved using dissipation-free spin currents.In this way, vortexbased EPR microscopy could be implemented without any externally applied magnetic field.High spin sensitivity stem from the small mode volume of the gyrotropic resonances which is independent of material parameters such as the saturation magnetization (cf.Fig. 1i).This is an important advantage of the vortex probe sensor compared to saturated ferromagnetic nano-objects, encouraging the use of ultra-low damping ferrimagnets that usually come with reduced saturation magnetization.This property, together with the possibility of using high frequency vortex modes, makes it possible to reach very large spin-magnon couplings.For example, a R = 100 nm YIG disc would allow detecting single spins with typical relaxation times of γ s ∼ 10 MHz.
Our approach is also potentially very interesting to increase the weak interaction between superconducting microcircuits and spin qubits, e.g., molecular qudits based on single rare earth ions that offer transitions between tunnel-split ground state doublets with high spin 74 .Using low-damping magnetic vortices could make it possible to read the state of individual spin qudits, which was previously considered unattainable 61 .Finally, we have also shown how to numerically normalize any magnon mode in ferromagnets of arbitrary size and shape.This can be used to calculate the zero-point magnetization fluctuations in confined nanomagnets including homogeneous magnon modes but also spin textures like domain walls, vortices or skyrmions.

A. Micromagnetic simulations
We use the finite difference software MUMAX3 75 to solve the time-dependent Laudau-Lifshitz-Gilbert equation for a given sample geometry and material.The relevant material parameters are the saturation magnetization M sat , the exchange stiffness A and the Gilbert damp-ing α (CoFe: M sat = 1.9 × 10 6 A/m, A = 2.6 × 10 −11 J/m and α = 10 −3 .Py: M sat = 0.86 × 10 6 A/m, A = 1.3×10 −11 J/m and α = 1×10 −2 .YIG: M sat = 0.14×10 6 A/m, A = 0.37 × 10 −11 J/m and α = 10 −5 ).In general, a disc of radius R and thickness t (or sphere of radius R) is simulated within a box with section 2R×2R and thickness large enough to calculate the relevant stray fields at 100 nm above the surface of the ferromagnet.The lateral size of the cells (with volume v cell ) is kept below o approximately equal to λ ex = 2A/µ 0 M 2 sat .Magnon dynamics are simulated using a dc biasing field B ap applied along a given direction, e.g.x in Fig. 1d.B ap is essential to obtain ferromagnetic resonances in the case of the saturated sphere but is unnecessary when dealing with vortices that do also resonate at B ap = 0.A sinusoidal time-dependent perturbation field B = βsinc(ω cutoff t) is applied perpendicularly to B ap in the case of the sphere (ŷ direction) and parallel to the disc plane (x, y) in the case of vortices.This is equivalent to exciting all spin-waves at frequencies below ω cutoff .We identify the magnon modes by calculating the numerical FFT of the resulting time-dependent spatially-averaged magnetization (see, e.g., Fig. 2c).We can write the dynamics at each cell for every mode, either a Kittel (K) or a vortex (v) one, as with r n the cell position, A(r n ) the amplitude and ω ξ their frequency.
To calculate the spin-magnon coupling we simulate the stray field, B rms (r j ), generated by the (zero-pointfluctuations of the) vortex magnetization from which the coupling can be obtained: For this purpose, we calculate the magnetic response on resonance using a perturbation field B = βsin(ω ξ t).Here ω ξ is the mode frequency and β must be low enough to keep the system in the linear response regime, i.e., β < 200 -500 nT.By doing so, we obtain, on the one side, the time-dependent vector magnetization at each cell n inside the ferromagnet.The amplitudes A in (1) depend on the perturbation strength β.Therefore, we must normalize them to having the single-magnon amplitudes.Using the usual formalism for magnon quantization the magnetization vector is normalized as 76 Here, M z is the spatial-averaged component of the magnetization along B ap in the case of the sphere and out-ofplane in the case of vortices.A α (r n ) are the perpendicular (the in-plane for the vortex case) amplitudes in (1) while δ x (r n )−δ y (r n ) is the phase difference between these two components.On the other side, we obtain the stray magnetic field resulting at every spin position r j outside the ferromagnet B stray (r j ) = B dc stray (r j ) + B ac stray (β, t; r j ) where we have split the static and the time-dependent components and the dependencies on position, time and perturbation field β are highlighted.The total zero-point field fluctuations are given by B rms (r j ) = ΛB ac stray (r j ), which is independent of β.Finally, we can calculate the contribution of B rms able of inducing spin transitions.The latter is given by B rms (r j ) = B 2 rms,1 + B 2 rms,2 where B rms,1 and B rms,2 are the components of B rms perpendicular to B tot (r j ) = B dc stray (r j ) + B ap (see Section Results).Inserting B rms (r j ) into ( 2), the positiondependent spin-photon coupling can be obtained.We emphasize that Fig. 1f and i show the spatial dependence of B rms (r j ), i.e., the part responsible for inducing spin transitions only.
The resonance window is also calculated numerically.For this purpose we consider only those spins for which the energetic criterion is satisfied, i.e., ω s = γ e B tot = ω ξ .This resonance window is enlarged by the broadening of both the magnonic resonance and the spins.In general, it will be reasonable to consider those spins that satisfy: with γ = max(γ v , γ s ).Finally, the average magnon-spin coupling per spin can be calculated as [Cf.Eq. ( 2)]: where the sum considers only those cells within the resonance window.This is to say, ( 5) is the sum of the position dependent coupling divided by the square root of the total number of spins that satisfy the resonance condition.

B. Electromagnetic simulations
In the case of the superconducting resonators, the spatial distribution of B rms,s is calculated using the finiteelement code 3D-MLSI 77 .This software solves London equations in a superconducting circuit with given dimensions and London penetration depth λ L .In the simulations we set λ L = 90 nm for Nb and a flowing current i rms .From here, 3D-MLSI allows calculating the spatial distribution of supercurrents and the resulting B rms,s (r j ).Under the particular geometry described in Fig. 1a, all spins satisfy the resonance condition and B rms,s (r j ) = |B rms,s (r j )|.The spatial-dependent coupling is given by Eq. (2).The experimental setup for spin detection discussed in the main text involves a LC resonator coupled to a transmission line.Simultaneously, the disc is coupled to the LC resonator (mounted on top of the inductor), as shown in Fig. 4 in the main text.The stabilized vortex scans the spins.The spin detection relies on the coupling of the latter to the vortex, resolved through a transmission experiment, that we model here.
An input coherent signal is sent through the transmission line, and the output is measured after interaction with the complete system: LC resonator + disc + spins.For convenience and clarity, we illustrate the coupling topology in Fig. 5.We observe that the LC is directly coupled to the TL.Thus, employing input-output theory, this transmission is given by: Here, ⟨r in,out ⟩ represents the averages of the input and output right-moving fields, κ quantifies the coupling to the TL, and χ a is the LC response, which can be computed as: Here, ⟨r in ⟩ = α in e −iωt .Hence, according to the inputoutput theory, the transmission reduces to the calculation of the LC response function in this case.To compute the LC response function, we employ a coupled mode model, depicted in Fig. 5, where the LC, vortex, and spin are modeled as harmonic modes coupled among them-selves.The dynamics are given by: Here, ⟨a⟩, ⟨a V ⟩, and ⟨a S ⟩ represent the modes for the LC resonator, vortex, and spins, respectively.The remaining parameters include the LC-vortex coupling, g LC , the (collective) spin-vortex coupling, G, and the line-LC coupling, κ.Finally, γ LC denotes the intrinsic losses of the resonator.The vortex and spin dissipation, γ v and γ s , are also taken into account in the model.Equation (8) represents a linear set of differential equations that can be solved by moving to the interaction picture ⟨a⟩ → ⟨a⟩e −iωt .This transformation yields an algebraic set of three coupled linear equations, which admit an analytical solution from which ⟨a⟩ can be obtained.Consequently, S 21 in Eqs. ( 6) and ( 7) is given by: The typical limiting cases of G = g LC = 0 (just the resonator) or the resonator plus vortex G = 0 are easily recognized.
The resonance frequency of the superconducting cavity is given by ω LC = 1/ √ LC with L and C the inductive and capacitive components of the circuit, respectively.L can be tuned by flux-biasing a superconducting quantum interference device (SQUID) coupled in series to the inductive part of the resonator 78,79 .In this way, ω LC can be adjusted so that ω LC = ω v0 .The latter is true even if a dc magnetic field is applied along the disc plane.Under these circumstances, the vortex resonance frequency ω v0 varies slightly (a few 10 MHz at maximum) as shown in Fig. 3 in the main text.On the other side, κ is fixed by design and typically takes values within a few kHz up to ∼ 100 MHz.The value of κ will determine the maximum quality factor of the resonator Q = ω LC /γ LC .In our simulations we use κ/2π = 30 MHz and Q = 10 3 which is feasible.
The magnetic disc behaves as a resonator with characteristic frequency ω v0 set by the aspect ratio and saturation magnetization M sat .Its loses are mainly given by the Gilbert damping α of the ferromagnetic material but they also depend on the particular excited mode.We determine both ω v0 and γ v by means of micromagnetic simulations with MUMAX3.Here, we use M sat , α and exchange stiffness A according to literature values for each material, as summarized in the Methods Section.
The coupling between the vortex and the superconducting resonator is calculated as described in Ref. 67,68 and briefly summarized here.We first calculate the zero point current fluctuations flowing though the LC resonator: with Z 0 the impedance.Next, we estimate the modulus of the zero point field fluctuations B rms,s produced by i rms at the vortex center position.B rms,s depends on the thickness and width of the superconducting line.Finally, the coupling is calculated as where V is the volume of the magnetic disc and ∆M is the amplitude of the volume averaged magnetization modulation when the vortex is excited with a varying magnetic field of amplitude b rms,s at frequency ω v0 .The validity of this equation has been demonstrated by comparing numerical simulations based on (11) with experimental coupling values from transmission and cavity measurements with Py nanomagnets 73 .
In our simulations we use 3D-MLSI to estimate the distribution of supercurrents in the superconducting circuits and the resulting magnetic fields.Here, we have assumed superconducting Nb circuits with thickness 50 nm and width 100 nm.The impedance of the LC resonator can be set by design and is assumed to be Z 0 = 10 Ω.By doing so, in the particular cases described in the main text we obtain g LC /2π = 6.8 MHz for the 400 nm ×60 nm Py disc and g LC /2π = 1.4 MHz for the 200 nm ×60 nm YIG disc.
Finally, we assume that each spin is resonant with the vortex at frequency ω v0 = ω S = γ e B tot with B tot = B ap + B stray the total static field at each spin position.G is calculated as described in the main text.

FIG. 1 .
FIG. 1.Comparison between electromagnetic cavities (coplanar waveguide) and magnonic resonators (saturated sphere and magnetic vortex).(a) (d) and (g): Scheme of the resonators.Bap is the external bias field (big blue arrow).Btot = Bap +Bstray is the total field (black arrows), with Bstray the demagnetizing field and Brms the zero point fluctuation field (gray arrows inside the red dashed circle).The yellow regions are the resonance windows.Below we plot the spatial distribution of Brms created by zero point current fluctuations in the electromagnetic cavity (b) and by the zero point magnetization fluctuations in a Py sphere (e) and a Py disc (h) of similar volume.Bottom panels (c, f and i) are the calculated Brms at 3 nm above the resonators at the position of the white dot in panels b, e and h.Simulations evidence the expected dependence Brms ∝ V 1/2 .

FIG. 2 .
FIG. 2. (a) radius of the vortex core rv vs t for R = 200 nm discs of different materials.rv is nearly independent of R. (b) Spatial distribution of the resonance frequency for free S = 1/2 spins for a YIG (top) and CoFe (bottom) discs with dimensions indicated.The latter is calculated as ωs = γeB demag .(c) Frequency spectrum of the same YIG (top) and CoFe (bottom) discs.The insets are the spatial FFT of the out-of-plane time-dependent magnetization of the different modes.

B 3 )FIG. 3 .
FIG. 3. (a) Normalized displacement of the vortex core d upon an external in-plane bias field B app∥ and (inset) spatial distribution of the magnetization of a R = 50 nm Py disc under B app∥ = 100 mT.(b) Resonance frequency for the v0 mode vs the external in-plane B app∥ (left) and out-of-plane field B app⊥ (right).In all panels, simulations are made for t = 60 nm and different radius (in nm) and materials as indicated in the legend.

2 FIG. 4 .
FIG. 4. (a) ⟨g⟩ vs t for different materials and radius (nm) indicated in the legend.(b) A Py disc is located over a superconducting LC resonator coupled to a transmission line (TL) and the S21 parameter is measured.The density plot shows the difference between the transmission at resonance with and without a small drop of DPPH at d = 52 nm from the disc center.(c) Calculated S21 with no disc, with a bare Py disc and for a disc containing a drop of ∼ 2000 DPPH spins.The inset shows S 21,bare − S21,DPPH in dB with the horizontal axis and scale being the same as in the main panel.(d) Transmission for a YIG disc (under identical conditions as in panel a) with B app⊥ = 5 mT.S21 is calculated with no disc, with a bare YIG disc and with a disc containing a single spin S = 1/2 at 3 nm over the surface with linewidth γs = 10 MHz and 1 Mz.The inset shows S 21,bare − S 21,S=1/2 in dB with the horizontal axis and scale being the same as in the main panel.

FIG. 5 .
FIG.5.Sketch of a transmission experiment.The input signal is sent through a transmission line interacting with the LC-circuit.In this sketch the modes (the LC-resonator, the vortex and the spins) are represented as circles, the couplings and dissipation channels are indicated.The dynamics of the model is given in Eq. (8).