Nanometre-scale probing of spin waves using single electron spins

Correlated-electron systems support a wealth of magnetic excitations, ranging from conventional spin waves to exotic fractional excitations in low-dimensional or geometrically-frustrated spin systems. Probing such excitations on nanometre length scales is essential for unravelling the underlying physics and developing new spintronic nanodevices. However, no established technique provides real-space, few-nanometre-scale probing of correlated-electron magnetic excitations under ambient conditions. Here we present a solution to this problem using magnetometry based on single nitrogen-vacancy (NV) centres in diamond. We focus on spin-wave excitations in a ferromagnetic microdisc, and demonstrate local, quantitative, and phase-sensitive detection of the spin-wave magnetic field at ~50 nm from the disc. We map the magnetic-field dependence of spin-wave excitations by detecting the associated local reduction in the disc's longitudinal magnetization. In addition, we characterize the spin-noise spectrum by NV-spin relaxometry, finding excellent agreement with a general analytical description of the stray fields produced by spin-spin correlations in a 2D magnetic system. These complementary measurement modalities pave the way towards imaging the local excitations of systems such as ferromagnets and antiferromagnets, skyrmions, atomically assembled quantum magnets, and spin ice.

. a, Reference frame xyz for the disc and x y z for the stray-field probe.
b, Field-dependence of the g(ω α,β ) rate for various magnetic fields, following Supplementary Eq.
(33) and using the parameters extracted from the fit in Fig.4b

Supplementary Note 1. Setup and sample
Our experiments were performed on a type IIa diamond grown by chemical vapor deposition (by the company Element 6) measuring 2x2x0.05mm 3 . We studied NV centres formed by N 15 ion implantation at 18keV and a density of 30/µm 2 and subsequent annealing for 2 hours at 800 • C, yielding NV centres at an estimated ∼50 nm below the diamond surface.
The magnetic fields used to control the NV centre spin state and to drive spin-wave excitations in the disc were generated by microwave (MW) currents. These currents were oil-immersion objective, resulting in the geometry depicted in Fig. 1a of the main text.
Laser pulses for optical spin initialization and readout were generated by an acousto-optical modulator (AOM) in double-pass configuration. The NV centre spin state was read out by integrating the first 550 ns of photoluminescence during a laser pulse. To avoid melting of the Permalloy disc, we limited the laser power to 250 µW.
The MW bursts used to control the NV spin state and to excite spin-waves in the disc were generated by modulating the output of a MW source (Rohde&Schwarz SMB100a or Agilent N5183A) using pulses generated by an arbitrary waveform generator (Tektronix AWG520) input into a MW switch (Minicircuits ZASWA-2-50DR+ or RFLambda RFSPSTA0208G).
To control the NV centre along two orthogonal axes (see Fig. 3a  For all the experiments, a static magnetic field was applied that separates the two NV centre spin transitions, which allowed the individual addressing of a target transition by MW pulses at the associated ESR frequency. To assure good optical spin contrast, the static field was carefully aligned with the NV centre crystal axis, using a procedure described in section Supplementary Note 7.

Supplementary Note 2. Numerical simulations
In this section we describe the numerical calculations of the spatial magnetization profile used in this work. Fig. 1f of the main text presented such calculations for the static magnetization of the disc. Fig. 2b of the main text shows the numerically calculated ferromagnetic resonance of the disc, and in Fig. 3 of the main text we presented numerical calculations of the time-averaged decrease in longitudinal magnetization of the disc.

A. Static magnetization and conventions
Micromagnetic simulations have been performed with the open source software OOMMF [16] running on the Harvard Odyssey cluster. For all the results presented in this work, the magnetic properties of the Permalloy disc have been simulated imposing a spatial discretization of 5×5×30 nm 3 . As in previous works, [20,21] we chose a saturation magnetization for Permalloy of M s = 800 kA/m, an exchange coupling of A=10 −11 Jm −1 , a Gilbert gyromagnetic ratio γ = 2.21 × 10 5 m/As, a damping constant α = 0.005. For the calculation of the static magnetic properties, the value of α has been increased to α = 0.95 to favour convergence. [16] Since the magnetic properties of a microdot are hysteretic [20], all measurements were done by first applying a large field and then sweeping the field down. To obtain Fig For the calculation of the stray field created by a certain magnetization pattern (see e.g. noting that even having perfect knowledge of the stray field B(r 0 ) in the entire plane at fixed distance d from the magnetic film, a unique reconstruction of the static spin structure is not possible, since any additional spin texture S(ρ) for which ∇ · S(ρ) = 0 will not affect the field. Accordingly, the kernel matrix present in eq. (5) of the main text is not invertible.

B. Uniform dynamics
Here we describe how we calculate and define the ferromagnetic resonance (FMR) frequency of the disc as a function of the external magnetic field B ext . We computed the linear, frequency-dependent, spatially-averaged, transverse magnetic susceptibility [15] χ ⊥ (ω) by applying a spatially uniform magnetic pulse h(t) of the duration ∆t = 50 × 10 −12 s oriented orthogonally to the plane of the disc (resembling the direction of the field generated by our coplanar waveguide). The following time-evolution m ⊥ (t) of the spatially integrated transverse disc's magnetization (transverse to the applied field and in the disc's plane) was then recorded for a total time of 20 ns, at 5 × 10 −12 s time intervals. We computed the field-dependent χ ⊥ (ω) via a Fast Fourier transform F , using the algorithm: In formulas, for the time evolution of the local (site-index i) magnetization m i (t), we assumed the following ansatz: Here, u ,i is the unit vector locally (index i) parallel to the equilibrium magnetization, and m eq ,i is the value of the local equilibrium magnetization. Furthermore, η labels the two axes orthogonal to u ,i and δm ,i (t) is the non-linear change of the longitudinal magnetization.
Provided with the previous assumption, we extracted from the numerical results the quantity δm ,i (t) by using the algorithm: The quantity δm ,i (t) is then time-averaged to compute the stray field along the NV axis ∆B (ω), to which our measurements are sensitive. The quantity ∆B (ω) is then casted into In this section we describe the fitting procedure for the Fano resonances observed in Fig. 2 of the main text. We use a phenomenological model by considering the interference between the response of a linear system and its driving force. In particular, we use the analogy between the response of a ferromagnet to a driving field and the response of a simple mechanical harmonic oscillator. In both cases energy is stored at resonance; the coherent precession of the two spin components transverse to the equilibrium axis of the ferromagnet is analogous to the periodic transformation of kinetic into potential energy. Furthermore, energy dissipation occurs via spin damping in the ferromagnet and through friction in the oscillator. Finally, both systems are characterized by a phase delay in their response upon driving. Let's consider the equation of motion for a simple mechanical oscillator: where m is the particle's mass, F (t) the driving force, β the friction, ∆ the resonance frequency and x(t) the time-dependent coordinate. The oscillator's response in frequency space reads as: where χ(ω) is the dynamical susceptibility. Similarly, the transverse dynamical magnetic uniform susceptibility χ ⊥ (ω) for a ferromagnet with short-range exchange interactions reads: [17] where W is the FMR line width, ∆ the FMR frequency and g(S) a constant prefactor which depends on the spin quantum number. A Fano interference is created in the mechanical system described by (5) when considering the total response: which implies a quadrature summation of the total normalized output's amplitude: x tot (ω) F (ω) = 1 + r 2 (ω) + 2r(ω) cos(θ(ω)).
Supplementary Eq. (8) is the mechanical equivalent of the expression we have used to model the enhanced normalized Rabi frequency. In particular, in the magnetic case, by calling the total transverse field at the NV i-site b i (ω) and the external driving field b D (ω) we obtain: where φ i is an additional frequency-independent phase factor motivated by the fact that the stray field created by the ferromagnet varies according to the different NV i-site. The terms θ(ω) and r(ω) can be obtained by taking argument and norm of the complex susceptibility in Supplementary Eq. (6). It is important to note that in order to model our field-dependent normalized Rabi curves (main-text Fig. 2b), one has to impose for the FMR resonance the Kittel-like expression ∆ = γ B ext (B ext + A) and for the frequency where D is the zero field splitting and B is the projection along the NV axis of the ferromagnet's static stray field (as measured in main-text Fig. 1c).
From the fits shown in Fig. 2b of the main text, we extract W/h = 0.2 GHz.
Supplementary Note 4. AC Stark effect in off-resonant detection scheme (Fig. 3 of main text) In this section, we analyze the phase shift imparted on the NV spin state by the AC Stark effect caused by the off-resonant driving field in the measurement scheme of Fig. 3A of the main text. The goal is to estimate if the Stark shift contributes significantly to the measured B ef f in Fig. 3 of the main text, taking into account that the AC drive field can be enhanced by the AC field generated by the ferromagnet as was shown in main-text Fig. 2.
The Stark shift, and the related Bloch Siegert shift, [5,6] is usually treated within an effective two-level model. [4,5] However, in our experiments, we prepare the S = 1 spin of the NV centre in a superposition of m s = 0, −1 states, while the off-resonant driving has a frequency in the vicinity of the m s = 0 ↔ +1 transition. Therefore, the three-level nature of the system has to be retained.
To discuss the effective field which is picked up by the NV centre as a result of the Stark shift, we employ the time-dependent Schrieffer-Wolff formalism. [7] The method outlined below can be readily generalized to other multilevel quantum systems.

A. Time-Dependent Schrieffer-Wolff formalism
The time-dependent Schrieffer-Wolff formalism is discussed in detail in Ref. 7. Here we briefly recall the main results. We assume an Hamiltonian composed of two parts, a timeindependent diagonal part H 0 and a time-dependent non-diagonal part H nd (t): In order to derive an effective Hamiltonian H ef f which retains up to second order the effects of the perturbation H nd (t), we change the quantum basis by applying a unitary In the Schrieffer-Wolff formalism, a series expansion ofŜ (t) leads to the expression for H ef f , which reads: [7] The matrix S(t) can be computed from the condition that eliminates the first-order terms O(H nd ), namely: [7] H nd (t) + Ŝ (t), H 0 + iŜ(t) = 0.
We will now use this formalism to obtain an expression for the time-independent part of H ef f for an NV centre in an off-resonant driving field.
where S = 1 spin matrices are used, and we have defined h = γB z and h 2 = γB 2 . We assume forŜ(t) the following 3x3 matrix representation, which satisfies theŜ(t) = −Ŝ(t) † condition:Ŝ Provided with Supplementary Eq. (12), the previous expression forŜ(t) defines a set of two linear non-homogeneous differential equations: where we have defined ω ± = D ± h, which are the spin transition frequencies in the absence of the Stark shift. We compute the solution to Supplementary Eq. (15) imposing S 1,2 (t = 0) = 0, i. e.,Û (t = 0) =Î. We obtain: The resulting time-independent partH ef f of the Hamiltonian (17) has the following form: Note that this effective Hamiltonian contains, in second order, off-diagonal terms that couple directly the m s = ±1 subspace. However, in our experiments, the degeneracy of the m s = ±1 states is lifted by the static external field h, which represents the dominant energy scale of the subspace. For this reason, as long as h h 2 , the spin dynamics will mainly be governed by the diagonal elements ofH ef f . In the limit h h 2 we are therefore allowed to write: which describes the shift of the energy levels of the NV spin in the presence of an off-resonant driving field.

C. Analysis of the Stark effect in the off-resonant detection scheme
We now consider the magnetometry sequence shown in Fig. 3a of the main text. The first π/2 pulse creates a superposition of m s = 0 and m s = −1 states, and all further NV-pulses are also given on the m s = 0 ↔ −1 transition. The normalized contrast is given by: where U x,y (π/2) denotes the operator for a π/2-pulse around the x, y axis. We obtain the following result: For illustrative purposes, we now consider two limiting cases of small detuning. If we apply off-resonant driving near the m s = 0 ↔ +1 transition at ω = ω + + δ, similar to Fig. 3b-c of the main text, we get for small δ.
On the other hand, if we would have applied driving near the m s = 0 ↔ −1 transition at which is the same as obtained from a two-level treatment of the NV-centre [5]. In this case, the Stark shift is twice that of Supplementary Eq. (21) because the driving is applied near the transition frequency of the two states forming the superposition. In this section, we estimate the magnitude of the Stark effect at the frequency of the FMR in Fig. 3 of the main text. The Stark effect quickly diminishes with increasing detuning, which is visible in Fig. 3c of the main text in the frequency range just above the 0 ↔ +1 transition over the entire magnetic field range. However, as we observed that the drive field may be enhanced by the spin-wave field (Fig. 2   To see if the Stark shift corresponds to a positive or negative effective field B ef f , we need to consider the Hamiltonian in Supplementary Eq. (19). As an example, we assume as in our experiments that ω = ω + + δ, with δ > 0, which leads to A < 0 while B ≈ 0. Suppose the NV is prepared in a arbitrary superposition Ψ(t = 0) = a|0 + b| − 1 + c| + 1 . The time evolution in a field h and with A < 0 will read as: effective field if we create a superposition of the 0,-1 states (with a, b = 0, c = 0) because the field h has changed to h + |A|. This positive Stark field corresponds to the experimental case. On the other hand, the effective field is negative for a superposition a, c = 0, b = 0, because the field h has changed to h − 2|A|. Therefore, the observation of the Stark shift, which as mentioned before is visible in Fig. 3c of the main text in the frequency range just above the 0 ↔ +1 transition, allows us determine the sign of B ef f in both the measurements in Fig. 3b-c of the main text (these measurements used the exact same pulse sequence).
Supplementary Note 5. Normalization procedure for the off-resonant detection scheme ( Fig. 3

of main text)
In this section we describe the normalization procedure used to obtain Figs. 3b,d of the main text.
We use the measurement scheme shown in Fig. 3a of the main text, in which we apply the final π/2-pulse along the x-or y-axis and subsequently read out the spin-dependent PL (P x and P y resp.). During the same measurement, we also apply two normalization sequences which are the same as in Fig. 3a of the main text except that we turn the MW off and apply the final π/2-pulse along the x and −x axis. Supplementary Fig. 5a shows the raw data of these four measurements. The normalization sequences yield the minimum and maximum PL (P min and P max resp.) which we use to obtain the x and y spin expectation values ( Supplementary Fig. 5b) according to i = 2 P i −P min Pax−P min − 1, where i = x, y. From the expectation values, we calculate the phase φ of the superposition using φ = arctan( y x ) (Supplementary Fig. 5c). We express the final signal in terms of an effective field B ef f = φ/(γT ), and divide by the square of the driving field |b d | 2 to correct for the frequency-dependent delivery of microwaves through our setup (Supplementary Fig. 5d).
We independently characterized |b d | 2 by measuring the Rabi frequency of NV ref at constant MW-source power as a function of the ESR frequency, which we tune using B ext ( Supplementary Fig. 6). The inset of Supplementary Fig. 5c illustrates the effect of the frequency-dependent power on the measurement of φ: dips/peaks in φ clearly correspond to dips/peaks in the Rabi frequency, motivating the normalization by |b d | 2 .
To validate the procedure of normalizing B ef f by the square of the drive field |b d | 2 (Supplementary Fig. 5d), we studied the dependence of B ef f on the power of the MW source R, since |b d | 2 ∝ R. Supplementary Fig. 7 shows the linear scaling of the signal with MW power.

Supplementary Note 6. Stray-field characterization of magnetization and spin noise
In this section we first describe the fitting procedure used to extract the relaxation rates from the measurements in Fig. 4a of the main text. Then we describe the model linking the relaxation rates to the spin-noise created by the disc. The following two sections are complementary to what discussed in the Methods section of the main text.

A. Noise probed by an NV centre
To describe the NV-spin relaxation we use a rate-equation model [13]: where P(t) describes the populations of the three NV-spin states as a function of time. In our fitting procedure, we imposed W 1,−1 = 0. This is validated by noting that magneticfield noise does not directly couple the m s = −1, 1 levels. In addition, we observed that the relaxation of NV A is dominated by magnetic-field noise (it has a much faster relaxation than far-away NV ref which we confirmed in a separate measurement to be ∼1/ms). The relaxation dynamics is therefore described by only two parameters: W 1,0 and W −1,0 .
In second order perturbation theory, the relaxation parameters are given by: [14] W α,β = α =β where the |α , |β are the spin eigenstates of the NV centre and ω α,β is the energy difference between the levels (ω α,β = (ω α − ω β )). The time-dependent magnetic perturbation at the NV centre due to magnetic-field fluctuations can be written as: whereĪ η is the η-component of the spin operator of the I = 1 spin of the NV and γ = 2π·28 GHz/T. It is easy to show that, by defining: we have: Besides probing the spin-noise spectrum at different frequencies, the two relaxation channels are formally identical.
Combining eq. (8) of the main text and Supplementary Eq. (27), we reach the following general expression for the characteristic relaxation rate of an NV centre in the vicinity of a thin magnetic film: Finally, note that the T -(temperature) and B ext -dependent spin fluctuations S η,η (ω α,β , k) can be related to the more well-known spin susceptibility χ η ,η (ω α,β , k) with the fluctuationdissipation theorem: [15]  Accordingly, we will compute the magnetic noise from a system composed by an infinitely thin ferromagnetic layer having a magnetic dipole density Γ. This procedure leads to the calculated noise spectrum and associated NV-relaxation rates shown in Fig. 4c of the main text (and to the corresponding solid lines in Fig. 4b of the main text.) We will give an estimate of the characteristic relaxation time including at first only in-plane transverse spin fluctuations, assuming that at the wavelengths k ∼ 1/d the exchange interactions are dominant. Assume the plane is magnetized along the z-axis. The susceptibility for a 2d ferromagnet along one of the axes transverse to z is given by: [17] χ y,y (ω α,β , k) = Γ −1 S W where D is the spin stiffness [17], W the width of the FMR excitation, ∆ its energy and S is the value of the local spin. In the previous expression the k in the denominator is expressed in m −1 . The value of D can be taken from previous works [18], namely D = 370 meVÅ 2 . Because of the demagnetization energy cost for out-of-plane spin fluctuations, we can in addition safely assume that χ y,y χ x,x and include in the calculations only the contribution of χ y,y .
In the case of Supplementary Eq. (31) the susceptibility has no φ k -dependence. We can therefore integrate out the φ k variable in the integral (29). We define the following θ NVdependent prefactor: We can now solve the integral (29), which gives: where the function G m,n p,q a 1 ,...,ap b 1 ,...,bq z is the Meijer G-function and d is the distance of the NV centre from the ferromagnetic film.
The previous expression formally holds true in the case of an infinitely large and thin magnetic film and it represents the contribution to the magnetic relaxation rate due solely to transverse spin fluctuations of the ferromagnet. On the other hand, longitudinal spin fluctuations take into account thermal intra-band spin-wave transitions; [17] due to the much smaller scale of the Zeeman with respect to the thermal energy, their contribution is essentially field-independent. We include them by summing up the constant g to the expression (33).
For a comparison with the experimental data we have used W/h = 0.2 GHz. This value was extracted from the fits of the Fano interferences in Fig. 2 of the main text with the damped oscillator response in Supplementary Eq. (6). We also vary ∆(H) as the Kittel law measured in our experiments for a 6 µm large disc and put S = 1/2 since fcc Permalloy has on average one Bohr magneton per site (M s · a 3 /4 ≈ µ B ) and a Landé factor g L ≈ 2 [22].
Finally, we used M s = 8 · 10 5 A/m as the saturation magnetization. [12] A plot of g(ω α,β ) with these parameters and d = 35 nm is in Supplementary Fig. 8b.
Supplementary Eq. (33) represents the noise produced by a two-dimensional ferromagnetic plane at a distance d. To illustrate the dependence on d, we shows fits to the noise measurements in Fig. 9 (same measurements as in Fig. 4b of the main text), in which we fixed d while fitting the parameters t and the field-independent offset g . When also releasing d, we obtain the fit for the full field-dependence of the spectrum g(ω α,β ) shown in Fig. 4b of the main text with d = 35(5) nm, in reasonable agreement with the ∼ 50 nm estimated from the implantation energy.
We now discuss the values of the fit parameter t, which underestimates the 30 nm expected for our evaporated Permalloy disc as can be seen from Supplementary Fig. 9. We note that several effects may play a role. Importantly, in our approximation the magnetic film is assumed infinitely thin, thus neglecting the smaller noise produced by spins that are located at a larger distance than the distance d from the NV centre. Furthermore, corrections due to the discretization of the spin wave spectrum associated with the finite thickness of the disc may also change the expected spin noise. Finally, a smaller saturation magnetization, fabrication-related imperfections, and/or oxidation of the disc may also lead to smaller spin noise.
In summary, our model for the spin noise captures the main physics of the problem, introducing a formulation of NV magnetometry in momentum space, which will be of use in future investigations of condensed-matter systems.

Supplementary Note 7. Photoluminescence-based field alignment procedure
In our experiments, we use the spin of the NV centre as an optically interrogated magnetometer. To assure good optical spin contrast, it is essential to align the applied magnetic field B ext with the N-V crystal axis [1]. In our experiments, all field sweeps are conducted by translating a permanent magnet along computed space trajectories that keep B ext aligned with the NV axis. In this section, we describe the procedure we used to find these trajectories, based on combining a model of the field generated by our magnet with a model of the NV-centre photoluminescence described by Tetienne et al [1]. In section Supplementary Note 7 C we describe experimental tests of the quality of the alignment.

A. Calculation of the photoluminescence
To apply B ext , we used cylindrical NdFeB magnets of variable dimensions. We calculate their space-dependent magnetic field profile with the open source Radia package [2]. Furthermore, we calculate the PL of an NV centre in a field using the rate-equation model and transition rates from Ref. 1. Combining these models yields characteristic plots for the space-dependent PL ( Fig. 10 left panel). The white dashed line marks the trajectory corresponding to perfect alignment. It is clear that the PL is very sensitive to misalignment whenever the field is ≈ 500 G or ≈ 1000 G. These field values correspond to the level anti-crossing of the ground and excited NV states, respectively. [1]. A direct comparison with experimental data around ≈ 500 G (Fig. 10, right panel) illustrates the quality of the model.
A transformation between the laboratory frame and the model frame is essential for calculating the magnet position corresponding to a given target magnetic field. The procedure to obtain this transformation is described next.

B. Determining the magnet trajectory along which the field is well aligned
To calculate the magnet trajectory along which the field is well aligned with the NV axis, we use the strong PL dependence on the magnet position in the B ext ≥ 514 G field range, which can be easily identified by the nodal point in the PL (Fig. 10 right panel). In particular, we first translated the magnet along orthogonal directions lying within N different constant-Z planes, recording the coordinates of a set of y i points, {i = 1, ..., N }, where the PL was found to peak. The y i represent a set of magnet positions in the laboratory frame for which the field is aligned. For each of these y i points we measured B ext,i by ESR. Provided with B ext,i , we numerically compute the corresponding x i positions of the magnet in the model frame. The transformation relating the y i to the x i is: where R is assumed to be a pure rotation matrix and t a translation vector. We obtain the matrix R via Wahba's method [3], which provides an expression for R based on the minimization of the following cost function: The algorithms for minimizing L (R) and computing R are described in Ref. 3. We then obtain the vector t by: Provided with the matrix R and the vector t, the laboratory frame coordinates of the magnet can be calculated for any value of B ext,i via (34). The next section describes experimental tests of the quality of the field alignment by measuring the field dependence of the ESR resonances of NV ref .
C. Test of the field alignment along the numerically calculated magnet trajectory In order to evaluate the quality of the field alignment along the numerically calculated maganet trajectory, we use single-NV vector magnetometry [11]. We measure the ESR frequencies ω ± of the m s = 0 ↔ ±1 transitions to determine the magnitude of the external field along (B z ) and transverse to (B ⊥ ) the NV axis.
In particular, we consider the NV Hamiltonian: where D is the zero-field splitting, and γ = 2π · 2.8025 MHz/G is the NV gyromagnetic ratio.
We substitute Supplementary Eq. (39) into Supplementary Eq. (38) and solve for λ 0 , B z , B ⊥ , using the measured values for ω ± . We obtain: The previous expressions rely on the value of D, which we extracted from the fieldindependent average D = (ω + + ω − )/2 to be D = 2π · 2.8707(1) GHz. In Fig. 11a we plot the extracted misalignment angle, defined as arctan(B ⊥ /B z ), for both NV ref and NV a .
A θ = 0 value for NV ref is solely due to a misaligned field of the permanent magnet, which is limited to the small value of θ ≤ 3 • (Fig. 11a). For NV a , the stray field from the Py disc causes the low-field increase of θ.

D. Determining the disc field from the ESR traces
The field created by the permanent magnet varies slightly between the different NVs studied in this work. In addition, the associated field gradient varies also with the magnitude of the external static field applied along the NV centre axis. By measuring the externally applied field at reference NV sites distant from the Permalloy disc, we measure a field gradient of  − ω NVa − )/γ. In the main text, we show the field dependence of B // , which includes the correction for the field-gradient.