Opportunities for long-range magnon-mediated entanglement of spin qubits via on- and off-resonant coupling

The ability to manipulate entanglement between multiple spatially-separated qubits is essential for quantum information processing. Although nitrogen-vacancy (NV) centers in diamond provide a promising qubit platform, developing scalable two-qubit gates remains a well-known challenge. To this end, magnon-mediated entanglement proposals have attracted attention due to their long-range spin-coherent propagation. Optimal device geometries and gate protocols of such schemes, however, have yet to be determined. Here we predict strong long-distance ($>\mu$m) NV-NV coupling via magnon modes with cooperativities exceeding unity in ferromagnetic bar and waveguide structures. Moreover, we explore and compare on-resonant transduction and off-resonant virtual-magnon exchange protocols, and discuss their suitability for generating or manipulating entangled states at low temperatures ($T\lesssim 150$ mK) under realistic experimental conditions. This work will guide future experiments that aim to entangle spin qubits in solids with magnon excitations.

Recently, several potential solutions to this challenge have been proposed by making use of boson modes as an information mediator. While photon-mediated NV-NV entanglement has been experimentally demonstrated over a meter and a kilometer length scales [10-12, 18, 20], based on indistinguishable single photon detection, its extension to two-qubit gates is still challenging due to its slow entangling rate as a result of its low success probability. It has been proposed, however, that the long-distance two-qubit gates can be realized by harnessing such entangled NV-center pair generation under both single-shot readout and local gates based on the measurement outcome [28]. This is possible if NV centers have access to quantum memories in the decoherence-free subspace [29], which survive during the multiple entangling attempts of NV centers that cause decoherence [13,[18][19][20]28]. Al-ternatively, as a means for extending NV-NV interaction on a wafer without needing single boson detection and with faster gate operations, hybrid quantum systems have been extensively studied where NV centers interface other bosonic systems [30][31][32][33]. In a carbon-nanotube-NV-center hybrid system [31], for example, it has been proposed to couple NV centers and phonon modes in a suspended carbon nanotube by injecting an electric current through the nanotube.
Here we present a practical and realistic hybrid quantum system to engineer NV-NV entanglement over micron length scales via on-and off-resonant magnon excitations at low temperatures (T 150 mK). The entanglement protocol in this hybrid quantum system is based on the strong coupling of NV spins to the magnon modes in yttrium-iron-garnet (YIG) nanodevices. Under a realistic geometry and accurately taking into account both dipole and exchange interactions, we obtain strong NVmagnon interactions and high entangling gate to decoherence ratio (GDR) in both an infinitely long YIG waveguide and a finite length YIG bar structure (see Fig. 1). Especially for the latter, we obtain NV-magnon cooperativity C 10 4 for on-resonance conditions and NV-NV GDR ≈ 10 3 under off-resonant magnon excitations for two NV centers separated by more than 2 µm. This leads to a usefully-fast entangling gate (relative to the qubit decoherence rate) at optically resolvable NV-NV separations. These values of GDR greatly exceed fidelities that were sufficient to demonstrate error correction on other platforms [53]. All of our results are obtained within a Hamiltonian formalism [54,55], which allows for semi-analytical expressions for the coupling in terms of the relevant experimental and geometrical quantities.
Finally, we explore and compare the calculated entanglement quality of both on-resonant transduction and offresonant virtual-magnon exchange entangling gate protocols, which we regard as another major focus in this work. We achieve this comparison by means of a numerical simulation of the Lindblad master equation taking into account two NV centers and a magnon mode near the resonance condition at finite temperatures. More specifically, we analyze and compare the entanglement negativity, fidelity, and degree of the Bell inequality violation for both cases under different parameters of the NV-magnon hybrid system. Notably, our results show that although the off-resonant protocols are robust at temperatures up to T ≈ 150 mK due to the absence of magnon occupation decay, the transduction protocol outperforms it due to its faster gate operations at lower temperatures if the magnon damping parameter is sufficiently small α (∆ω/ω µ )(1/4g µ T * 2 )[π/(π − 1)], with magnon frequency ω µ , NV center coherence time T * 2 , NVmagnon detuning frequency ∆ω, and NV-spin-magnonmode coupling g µ . Our calculations and analysis serve as a guide for future experiments to engineer on-chip longdistance entangling gates between NV centers mediated by magnons in ferromagnetic nanostructures.
In this article, we begin in Sec. II with the description of the Hamiltonian formalism for the dipole-exchange magnons coupled to NV centers. In Sec. III we calculate the full magnonic properties of a YIG waveguide interacting with NV centers. We obtain the NV-NV coupling strength, the entanglement rate, and the gate to decoherence ratio under the off-resonant NV-magnon interaction condition. Similarly, in Sec. IV we first calculate the magnonic properties of a finite length YIG bar. Secondly, we evaluate both NV-magnon on-resonant coupling strength and its cooperativity as well as the NV-NV coupling strength under the off-resonant condition. We provide for the latter the entanglement rate and the gate to decoherence ratio. Finally, in Sec. V we present a complete comparison between the transduction and virtualmagnon-exchange protocols in detail under different system parameters and physical conditions.

II. HAMILTONIAN FORMALISM OF DIPOLE-EXCHANGE MAGNONS AND NV-MAGNON INTERACTION
Here we outline the Hamiltonian formalism of dipoleexchange magnons coupled to NV centers providing a complete and accurate treatment of both magnetic dipole and quantum exchange interactions between the spins in YIG waveguides and bars with finite cross section. This is crucial in our study as the NV centers have eigenfrequencies typically on the order of gigahertz, thus interacting with the so-called dipole-exchange magnons in ferromagnets [50]; using simpler, less accurate magnon dispersion relations as in Ref. [34] leads to a substantial overestimation of the NV-magnon coupling. As illustrated in Fig. 1, we consider hybrid quantum devices where NV centers are placed on top of the YIG structures. Whereas multiple NV centers can be placed on top of the infinitely long YIG waveguide in a scalable fashion as shown in Fig. 1, in the following calculations we only focus on coupling two NV centers. The total Hamiltonian of our hybrid system is written as H = H NV + H m + H int , where H NV is the NV Hamiltonian, H m is the magnon Hamiltonian, and H int is the interaction Hamiltonian, H m = −µ 0 drH ext · M(r) + µ 0 2 drα ex (r)∇M : ∇M Here, D NV = 2π × 2.877 GHz is the zero-field splitting of the NV center,n NV is the unit vector along the NV main symmetry axis, S NV i is the spin-1 operator of the NV center labeled by i ∈ {1, 2}, γ = 2π × 28 MHz/mT is the absolute value of the electron gyromagnetic ratio, µ 0 is the vacuum permeability, H ext is the external magnetic field, M(r) is the magnetization with the constraint |M(r)| = M s (r) = M s F (r), M s = 245.8 mT/µ 0 is the YIG saturation magnetization, F (r) = 1 (0) inside (outside) the ferromagnetic structure, α ex (r) = α ex F (r), α ex = λ 2 ex = D ex /γµ 0 M s is the exchange-length squared, D ex = 5.39 × 10 −2 γ mT µm 2 is the YIG exchange constant, the double-dot product is defined as ∇M : ∇M = ∂ a M b ∂ a M b , r i is the position of NV i , G(r − r ′ ) = 1/4π|r − r ′ | is the Green's function, and we set = 1. We note that the first term in Eq. (2) is the Zeeman energy, the second term is the exchange energy, and the third term is the magnetic dipole energy. Inclusion of both the second and the third term in Eq. (2) results in the dipole-exchange magnons in ferromagnets.

III. INFINITELY LONG FERROMAGNETIC WAVEGUIDE
Here we consider the case of an infinitely long YIG waveguide with thickness, width, and length given by d, w, and l(→ ∞), respectively. The external magnetic field is applied along the YIG waveguide, H ext = H extẑ , and NV centers are positioned at height h from its top surface [see illustration in Fig. 2(a)]. The equilibrium magnetization is M 0 (r) = M sẑ F (r), for which its contribution in the interaction Hamiltonian Eq. (3) vanishes. The NV main symmetry axis is set to be parallel to the external magnetic field,n NV =ẑ, for geometrical simplicity. We further define the deviation from the equilib- ŷ is a small two-dimensional magnetization deviation. The linearized magnetization dynamics [56] are governed by the Hamiltonian equation of motion for m − (r) = [2γM s (r)] 1/2 a(r) and m + (r) = [2γM s (r)] 1/2 a * (r) using the magnon Hamiltonian H m up to quadratic order in the complex canonical variables a(r) and a * (r), where we have performed the Holstein-Primakoff approximation [51] and m ± (r) = m x (r) ± im y (r).
To obtain the normal magnon mode frequencies and the dynamical fringe field spatial profiles, we diagonalize the magnon Hamiltonian Eq. (2) by expanding the complex canonical variables assuming totally unpinned surface spins, i.e., Here, the basis functions are where , F Y (y) = Θ(y)Θ(w−y), and Θ is the Heaviside step function. As we consider the case where both the thickness and the width of the YIG waveguide are small, we restrict our discussion to the magnon mode subspace with (n, m) = (0, 0), which presents uniform magnetization deviations across the x-y plane and gives the lowest energy magnon band in the dispersion relation.
After writing H m up to the quadratic order in the complex canonical variables, applying the Bogoliubov transformation, and promoting the complex canonical variables to the quantum creation and annihilation operators, we obtain the diagonalized Hamiltonian (see Appendix B1) where ω k,(0,0) is the magnon energy and β k,(0,0) (β † k,(0,0) ) is the magnon annihilation (creation) operator satisfying . The coupling strength between magnon modes and NV centers can be obtained by applying the same Bogoliubov transformation in the interaction Hamiltonian Eq. (3). As we focus on external magnetic field values γH ext < D NV , the NV center's ground state and the first excited state are |g = |S z NV = 0 and |e = |S z NV = −1 , respectively. Up to the linear order in magnon creation and annihilation operators and using the rotating wave approximation (|ω k,(0,0) − ω NV | ≪ ω k,(0,0) + ω NV ), we obtain the interaction Hamiltonian (see Appendix B2) in the NV centers' subspaces spanned by is the dimensionless coupling between the NV center spin and the k-magnon mode, ρ i is the NV i 's position in the x-y plane, σ + NVi = |e i g|, and σ − NVi = (σ + NVi ) † . The virtual-magnon-mediated NV-NV interaction can be obtained via the Schrieffer-Wolff transformation [57] as where g eff is the effective NV-NV coupling strength, ω NV = D NV − γH ext is the transition frequency of |g ↔ |e , and we write g(k) = g(ρ i , k) assuming ρ 1 = ρ 2 . The above expression is valid when We note that this effective coupling strength g eff for the offresonant configuration does not depend on the temperature, as it is independent of the initial magnon number state |n m (i.e. from second order perturbation theory) even though the NV-magnon coupling strength matrix element is proportional to √ n m + 1 (see Appendix B4). In Fig. 2(b) we plot the NV center's transition frequencies and magnon mode frequencies as a function of the external magnetic field H ext , where we have assumed (d, w) = (20 nm, 120 nm) for the waveguide dimensions [58]. As we take the limit where the length of  9)] as a function of the NV-NV distance under ∆f = 3 MHz and ∆f = 10 MHz. The gray curve shows the coupling due to the direct magnetic dipole-dipole interaction between NV centers. The entanglement rate and the gate to decoherence ratio are shown on the right axis for T * 2 = 1 ms. Inset shows the time τ evolution of the entanglement negativity at T = 0 from the initial state |g 1|e 2 scaled by the Bell state negativity NB.
the YIG waveguide is infinity (l → ∞), the magnon mode frequencies form a continuum with its minimum denoted as ω min . At field H ext = H c , the NV center's lower transition frequency ω NV is detuned from the magnon dispersion minimum ω min by ∆ω = ω min − ω NV = 2π∆f = 2π × 3 MHz. Figure 2(c) shows the magnon dispersion relation near ω min and the wavenumber dependence of the dimensionless coupling strength g(k) at H ext = H c , ρ i = (d + h)x + wŷ, and h = 25 nm [see the cross marker in Fig. 2(d)]. The coupling strength also depends on the spatial position of the NV center relative to the YIG waveguide, which is shown in Fig. 2(d). As the dynamical fringe magnetic field generated by a single magnon is confined near the YIG device, the coupling strength is larger if the NV center is positioned near the YIG waveguide.
Under the off-resonant condition shown in Fig. 2(c), the NV centers on top of the YIG waveguide interact to each other via the exchange of virtual magnons. In Fig. 2(e), we plot the effective NV-NV coupling strength g eff [Eq. (9)] as a function of the NV-NV distance δz = |z 1 − z 2 | for both ∆f = 3 MHz and ∆f = 10 MHz cases represented by the red and blue dots, respectively. The coupling decays rapidly with detuning, which allows the entangling interaction to be switched off by increasing the external magnetic field from H ext = H c by ≈ 0.1 mT. We show that the calculated coupling strength is well explained by the analytical formula as shown by the solid red and blue curves in Fig. 2(e), where ξ 0 = D ex /∆ω is the spin correlation length and ωd = µ 0 γ 2 /(ξ 0 wd). The entangling gate rate ER = 4g eff /π and the gate to decoherence ratio GDR = 4g eff T * 2 /π are shown on the right axis, where a coherence time T * 2 = 1 ms of the NV center is used [7]. As we obtain GDR > 10 for 1 µm separated NV centers, we predict a useful and practical entangling gate.
To show that this system can manipulate the NV-NV entanglement, we perform a simulation using the Lindblad master equation. In the inset of Fig. 2(e) we plot the entanglement negativity [59] N at T = 0 as a function of the NV-NV interaction time after the preparation of the initial spin state in |g 1 |e 2 , where the negativity is normalized by the Bell state's negativity N B . As we obtain N > 0, we clearly demonstrate that the NV centers are entangled. If multiple NV centers are placed on top of the YIG waveguide (see Fig. 1), neighboring two-NV gates can thus be performed by locally changing the external magnetic field around the two NV centers to shift their transition frequencies relative to the minimum magnon mode frequency in the range ∆ω > 0. Alternatively, local electric field [60] or strain [61] can be used to shift NV centers' transition frequencies to avoid applying a local magnetic field at the underlying YIG location, the effect of which is discussed in Appendix K.
In Fig. 3 we plot the NV-NV entanglement rate and the gate to decoherence ratio as a function of the waveguide thickness d for different waveguide dimensions and NV centers' heights h. We assume a fixed NV-NV distance of 1 µm, (x i , y i ) = (d + h, w), and ∆ω = 2π × 3 MHz. The red (blue) solid curve shows the waveguide thickness d dependence of the ER and the GDR under the fixed aspect ratio w/d = 6 at h = 25 nm (5 nm), and the red (blue) dashed curve shows the dependence where the waveguide width is kept constant with w = 120 nm at h = 25 nm (5 nm). From these graphs we see that in order to make the entangling gate faster, one can either have the NV center closer to the YIG waveguide (diminishing h) or make the waveguide's cross-sectional area smaller. As for placing NV centers in proximity to the YIG waveguide, we note the common challenge of making high coherence NV centers near the diamond surface due to the surface noise known in the area of NV-based quantum sensing [62].

IV. FINITE LENGTH FERROMAGNETIC BAR
In this section we show that the NV-magnon coupling strength can be strongly enhanced under the magnon confinement effect of a finite length ferromagnetic bar. As the magnon mode frequencies are discretized for this case, the system allows us to control the NV levels to be on-and off-resonant to the magnon levels. Here, the interaction Hamiltonian Eq. (8) can be transformed into the form of the Jaynes-Cummings model [39,63], and the entangling gate schemes used in both quantum optics and circuit quantum electrodynamics can now be implemented in our hybrid quantum system [64][65][66].
We first obtain the NV-magnon interaction Hamiltonian for a finite length YIG bar using a similar procedure as done in Sec. III. For that, we first take the equilibrium magnetization to be M 0 = M s F (r)ẑ and approximate the x, y component of the resulting static demagnetization field in Eq. (2) to be negligible compared to its z component. Although there is also a finite static demagnetization field contribution in the interaction Hamiltonian Eq. (3), we verified that its value is small under the geometry parameters and NV center positions we consider.
Accordingly, we diagonalize the magnon Hamiltonian through the following expansion of the complex canonical variable where the z-directional basis function is κ Z p = pπ/l, and F Z (z) = Θ(z)Θ(l − z). As we consider the case with d, w ≪ l, we restrict our discussion to the magnon mode subspace with (n, m) = (0, 0). Considering z-directional modes with p = 0, 1, · · · , N , where p = N labels the highest z-directional wavenumber mode to be taken into account, and keeping terms up to the quadratic order in the complex canonical variables, we obtain a 2(N +1)×2(N +1) non-diagonal quadratic boson Hamiltonian. After applying the Bogoliubov transformation with the paraunitary matrix [54,56] and promoting the complex canonical variables to the quantum creation and annihilation operators, we obtain (see Appendix C1) In a similar way as in Sec. III, the NV-magnon interaction Hamiltonian can be mapped into the form of the Jaynes-Cummings model [39,63] (see Appendix C2) where g µ (r i ) ∝ √ ω M ω dwl [ω dwl = µ 0 γ 2 /(dwl)] is the coupling strength between the NV center spin and the µ-magnon mode in the unit of energy. As the magnon creation operator β † µ applied to the magnon number state |n µ gives rise to a factor of n µ + 1, we expect the on-resonant NV-magnon configuration to have n µ + 1 faster energy-transfer oscillations between the NV-center spin and the µ-magnon mode. However, at finite temperature, which can be thought of as a statistical mixture of different magnon-number states, these different-period oscillations will average out incoherently. Therefore, finite temperature does not improve the quality of NV-NV entanglement via magnon modes even though the mean magnon number n µ is larger, indicating that magnonmediated NV-NV entanglement needs to be performed at low temperatures T 150 mK (see Sec. V).
In Fig. 4(a) we plot the external magnetic field H ext dependence of the discretized magnon mode frequencies of a YIG bar with dimensions (d, w, l) = (5 nm, 30 nm, 3 µm). The neighboring magnon mode frequencies are separated from each other by over 2π × 10 MHz for modes with p ≥ 5, as shown in Fig. 4(b). At field H ext = H c , the NV center's transition frequency ω NV and the magnon mode frequency ω (005) are on-resonant. We plot in Fig. 4(c) the spatial distribution of the NV-magnon coupling strength g (005) at H ext = H c for a fixed NV center height h = 5 nm The white rectangle delimits the bar dimension, and the white cross mark represent the position of NV1 referred in (d). The corresponding cooperativity C (005) is shown on the right axis. (d) Effective NV-NV coupling strength g eff between two NV centers as a function of the NV-NV distance, where NV1 and NV2 are placed at r1 = (d + h)x + wŷ + (400 nm)ẑ and r2 = r1 + δzẑ, respectively. The red (blue) curve is calculated for (d, p) = (5 nm, 5) [(d, p) = (20 nm, 12)]. The entanglement gate rate (ER) and the gate to decoherence ratio (GDR) are shown on the right axis. In both cases aspect ratio is w/d = 6, length of the magnetic bar is l = 3 µm, and detuning is ∆f = ∆ω/2π = (ω (00p) − ωNV)/2π = 3 MHz.
[see Fig. 2(a)], and obtain g (005) ≈ 2π × 0.5 MHz depending on the NV center positions. With the Gilbert damping parameter of YIG α = 10 −5 [67] and the coherence time of NV centers T * 2 = 1 ms [7], we show on the right axis of Fig. 4(c) the corresponding single magnon µ-mode cooperativity [30,33] which is a dimensionless measure of the coupling. We emphasize that because this represents the single-magnonmode cooperativity, the temperature dependence only appears in α and T * 2 which for the purpose of our lowtemperatures analysis are assumed to be independent of temperature. We find C (005) 10 4 over a wide range of NV center positions, achieving the strong coupling regime for our hybrid quantum system. In contrast to Sec. III, where we have a translationally invariant infinitely long waveguide, here the position of the NV center along z-direction plays a major role in the coupling strength. Our calculations enable us to optimize both the coupling strength and the cooperativity in order to increase NV-NV entanglement efficiency in our system.
The virtual-magnon-mediated NV-NV interaction is calculated in a similar way as in Eq. (9) under the condition |g µ (r)| ≪ |ω µ − ω NV |, and we obtain with µ = (005) (see Appendix C3). In the same way as in the waveguide case, this virtual-magnonmediated coupling strength is independent of temperature. Here, the two NV centers are placed at Fig. 4(c)] and r 2 = r 1 + δzẑ, where δz is the NV-NV distance along the bar length. In Fig. 4(d) we plot g eff as a function of δz for the detuning ∆ω = ω (005) − ω NV = 2π × 3MHz, which could be produced by electric field [60], strain [61] or magnetic field deviation from H ext = H c . The corresponding entangling gate rate and the gate to decoherence ratio are shown on the right axis. Surprisingly, useful entangling gates for 2.2 µm separated NV centers with g eff = 2π × 90 kHz and GDR > 700 are predicted for this YIG bar system. This makes experiments more accessible in terms of the independent optical initialization and the readout of NV centers than the waveguide case.
We have also calculated these quantities for a less challenging to fabricate YIG geometry with (d, w, l, h) = (20 nm, 120 nm, 3 µm, 5 nm). The result is plotted as a blue curve in Fig. 4(d), for which we obtain GDR > 100 for the 2.2 µm separated NV centers. This result clarifies the significance of using the YIG bar structures to entangle two NV centers separated by a few micrometers. Moreover, the discretized magnon mode frequencies allows for controlling the NV center frequencies to be on-resonant to one of the magnon mode frequencies, which enables the entanglement of two NV centers via the transduction of energy quanta that we discuss in the next section. We also comment that it would be possible to control the NV-magnon coupling strength via parametric driving of the discretized magnon modes as studied in the cavity quantum electrodynamics [68] (see Appendix I).

V. TRANSDUCTION AND VIRTUAL-MAGNON EXCHANGE PROTOCOLS
In this section, we explore and compare two entangling gate protocols for our hybrid quantum system, onresonant transduction and off-resonant virtual-magnon exchange. Entanglement via the transduction protocol is simulated by controlling the NV center frequencies independently as illustrated in the left schematic of Fig. 5(a). For this case, the NV spins are initially prepared in the state |g 1 |e 2 , i.e., NV 1 (NV 2 ) is in its ground (excited) state. We first make ω NV 2 on-resonant to the µ-magnon mode frequency ω µ for a certain time τ var during which ω NV 1 is detuned from ω µ by δω =2π × 5 MHz. Second we swap the NV 1 spin state and the magnon state by making ω NV1 = ω µ for the swap gate time τ SWAP during which ω NV 2 is detuned from ω µ by δω. The total interaction time in this protocol is τ int = τ var + τ SWAP and is varied by changing τ var . The control of the NV centers' transition frequencies can be performed by applying a local magnetic field, electric field [60], or strain [61]. An alternative possibility of controlling the transition frequencies would be to use a periodic modulation of the external magnetic field [69,70] (see Appendix J). In contrast, in the virtual-magnon exchange protocol the NV centers' frequencies are both detuned from the µ-magnon mode frequency by ∆ω = ω µ − ω NV 1,2 = 2π × 3 MHz [see the right schematic of Fig. 5(a)]. After the preparation of the NV centers' spin state in |g 1 |e 2 , the whole system evolves over the interaction time τ int .
The time evolution of our hybrid quantum system for both protocols is simulated using the Lindblad master equation [30,71,72] at a finite temperature T considering two NV centers and a magnon mode µ, where is the Boltzmann constant, and ρ is the density operator. Here, the magnon damping parameter κ = αω µ is based on the dissipation term in the Landau-Lifshitz-Gilbert For the magnon mode contribution in the total Hamiltonian H(t), we only take into account the magnon mode with µ = (005), as this mode produces the dominant contribution in the NV-NV interaction as well as the magnon induced decoherence of NV centers in both protocols. As the NV center's longitudinal decay rate is much smaller than the transverse decoherence rate [6], we assume it to be zero in the simulation. The two NV centers are separated by 2.2 µm along the YIG bar length with r 1 = (d + h)x + wŷ + (400 nm)ẑ and r 2 = r 1 + (l − 800 nm)ẑ. We use the Gilbert damping parameter α = 10 −5 of YIG [67] and the NV center coherence time T * 2 = 1 ms [7]. In the upper two panels of Fig. 5(b), we plot the NV centers' excited state population p ie (i = 1, 2) and the magnon population n = n µ [µ = (005)] at the end of the transduction (on resonant) and the virtual-magnon exchange (off resonant) protocols as a function of the total system interaction time τ int at T = 70 mK. In the lower two panels we plot three different entanglement measures as a function of the interacting time τ int for each protocol. More specifically, we plot the entanglement negativity normalized by the Bell-state's negativity, the degree of the Bell inequality violation, and the fidelity to the target pure entangled states, which are given by the red, sky blue, and gray curves, respectively. The resulting states are entangled if N > 0, and one expects to observe the violation of the Clauser-Horne-Shimony-Holt (CHSH) form of Bell inequality if CHSH Violation > 0 [73,74] (see Appendix D1).
In Fig. 5(b) we first find that the transduction protocol is faster in gate operation as compared to the virtualmagnon exchange protocol. This is because the NVmagnon on-resonant coupling rate g µ ≈ 2π × 0.5 MHz is larger than the off-resonant NV-NV coupling rate g eff ≈ 2π×90 kHz. On the other hand, it is observed that the virtual-magnon exchange protocol results in larger amplitude oscillations in the NV centers' excited state populations and higher fidelity under the parameters and the temperature used in the simulation. This result is understood by a combination of two factors. First, the virtual-magnon exchange protocol only creates magnons virtually (with magnon population suppressed by g µ /∆ω due to the energy mismatch), thus being approximately insensitive to the magnon damping parameter. Secondly, the magnon damping rate αω µ is faster than the NV center's decoherence rate 1/T * 2 , and therefore there is more loss of information if a real magnon is excited. Nonetheless, in both protocols we predict entangled states can be manipulated and the violation of the Bell inequality will be observed.
To further compare the two entanglement protocols, we have performed simulations under multiple temperatures and have observed that the virtual-magnon-exchange protocol is more robust at higher temperatures up to ≈ 150 mK (see Appendix D2 and E). Moreover, we show that both protocols do not produce useful entanglement for T 150 mK due to the NV centers' dephasing from magnon number fluctuations of modes with µ = (005). We have also evaluated the decay contribution due to these magnon modes and have verified that this is negligible for temperatures T ≤ 150 mK for both upper and lower transitions of NV centers (see Appendix H). Interestingly, the transduction protocol improves more drastically at lower temperatures than the virtual-magnon exchange protocol. Based on the zero temperature analysis, we find an inequality for which the transduction protocol performs better (see Appendix D2) For the parameters used in this section, the transduction protocol is shown to outperform the virtual-magnon exchange protocol (with ∆ω = 2π ×3 MHz) if α 10 −7 . In Appendix D2 we provide phase diagrams in (α, 1/T * 2 )space for which protocol gives higher fidelity under multiple detuning values. Analytical expressions for the fidelity in the limit αω µ /g µ ≪ 1 and T * −1 2 /g µ ≪ 1 are also provided. To show that the magnon-mediated entanglement scheme can directly be extended to two-qubit entangling gates, we have also calculated an average gate fidelityF [75] as a square-root-of-iSWAP gate for the off-resonance protocol, and have obtainedF ≈ 0.88 at T = 70 mK (see Appendix F).
As for keeping the system at low temperatures T 150 mK, we note that the laser illumination and microwave irradiation on the system for the initialization, manipulation, and readout of NV centers may cause unwanted heating. Although YIG has been studied under microwave irradiations in superconducting qubit platforms [76] and color centers have been studied under laser illuminations in dilution refrigerator temperatures T < 100 mK [77][78][79][80], it would be important to minimize the average microwave irradiation and laser illumination power on the system to maintain the required low temperatures. Here, of particular interest is the possibility of cooling down the target magnon mode to its ground state in analogy to cavity optomechanics techniques [81][82][83][84][85], e.g., via the optomagnonic interaction [86] or via the coupling to NV centers [84,85]. For example in Fig. 5(b), we have observed that the mean magnon occupation number at the end of the on-resonant protocol is smaller than its thermal level [see n(τ int = 0) in the off-resonant protocol], which is reminiscent of the ground-state cooling of magnons and motivates future studies on the alternative cooling methods of the NV-magnon hybrid quantum system.
We also note that the small Gilbert damping parameter α = 10 −5 used in the current study may be optimistic for small YIG structures as the value is obtained from bulk YIG samples [67]. This is partially due to the nontrivial magnetic behavior at millikelvin temperatures of the gadolinium-gallium-garnet (GGG) substrates on which YIG is typically grown [87], which would be mitigated by employing a free-standing structure [88], and also due to the impurity relaxation mechanism in YIG [89]. However, with remarkable advances in recent magnonics research, it has been shown that the damping of thin YIG films can be improved considerably, e.g., with techniques based on a recrystallization of amorphous YIG into single crystals [90]. Additionally, we obtain a high cooperativity C ≈ 500 even with the larger Gilbert damping parameter α = 10 −3 as calculated from Fig. 4(c). We have further performed simulations with α = 10 −3 in Appendix G, and find that the entangled state can still be produced at T = 70 mK for the off-resonant protocol, although further optimization on the detuning frequency is needed to improve the quality of the entanglement in order to avoid the overlap of the NV centers' transition frequencies with the now broader linewidth of the magnon mode resonance (see Appendix G).

VI. CONCLUSION
We study hybrid quantum systems consisting of NV centers in diamond and magnons in ferromagnetic bar and waveguide structures. Based on the Hamiltonian formalism of the dipole-exchange magnons, we predict useful two-NV entangling gates over 1-2 µm NV-NV separations at finite temperatures. Transduction and virtualmagnon exchange protocols of entanglement are explored and compared under realistic experimental conditions. Although the transduction protocol is faster in gate operation, the virtual-magnon exchange protocol results in higher fidelity as the typical Gilbert damping parameter of YIG makes the magnons less coherent than the NV centers. We have obtained entangled state's fidelity F ≈ 0.81 for the transduction protocol and F ≈ 0.95 for the virtual-magnon exchange protocol at T = 70 mK. The virtual-magnon exchange protocol is also found to be robust against thermal magnon fluctuations, although the transduction protocol outperforms it close to zero temperature for αω µ T * 2 (∆ω/g µ )/[4(1 − 1/π)]. Calculations presented in this study help to implement optimal device geometries and entangling gate protocols in future experiments trying to entangle spatially separated NV centers using magnons in ferromagnets. The total Hamiltonian H of our hybrid system composed of NV centers and magnons is presented in Sec. II. We note that the interaction Hamiltonian H int can also be understood in terms of the dipolar tensorD(r − r ′ ): The magnetization dynamics governed by the Landau-Lifshitz-Gilbert (LLG) equation (without dissipation) is obtained by the Hamiltonian equation of motion with the following identification of the magnetization and the classical complex canonical variables following the Holstein-Primakoff transformation: where Here, a(r) and a * (r) are the complex canonical variables satisfying ∂ t a(r) = −iδH/δa * (r) and ∂ t a * (r) = +iδH/δa(r). The relation between a(r), a * (r) and M − (r), M + (r) is carefully chosen such that it satisfies the dissipationless LLG equation ∂ t M = −γµ 0 M × H eff . This is also consistent with the standard sign convention of the time evolution of the creation/annihilation operators a ↔â ∝ e −iωt and a * ↔â † ∝ e +iωt . Note that we use γ > 0 as the absolute value of the electron gyromagnetic ratio, so the electron gyromagnetic ratio is −γ. The Hamiltonian equation of motion gives and writing As the effective field is obtained by In the following discussions, we apply an external field along theẑ direction of Fig. 2(a), H ext = H extẑ , and for geometrical simplicity we take the NV main symmetry axis to be alongẑ axis, i.e.,n NV =ẑ. We will consider the case where the equilibrium magnetization is uniform across the ferromagnet and parallel toẑ, M(r) = M 0 (r) = M s (r)ẑ. Although in principle we need to obtain the equilibrium magnetization from the energy minimization of H m , in the infinitely long waveguide case M(r) = M s (r)ẑ holds as the field is applied along the direction where demagnetization factor is zero. In the finite length magnetic bar case, this is still approximately correct as in our setting the length l is much larger than both the width w and the thickness d. Here, m(r) = m x (r)x + m y (r)ŷ is a small twodimensional magnetization deviation. Now the deviation can be written as δM(r) ≈ m(r) − (m 2 (r)/2M s (r))ẑ.

Simplification of the magnon Hamiltonian
In the following calculation, we simplify the magnon Hamiltonian H m . We write where H Z is the Zeeman Hamiltonian, H ex is the exchange Hamiltonian, and H dip is the magnetic dipole Hamiltonian given by where the double-dot product is defined as ∇M : ∇M = ∂ a M b ∂ a M b . Firstly, we simplify the Zeeman Hamiltonian and the dipole Hamiltonian. Using M(r) = M 0 (r) + δM(r), we obtain Here, H dem is the demagnetization field Hamiltonian, H dip (2) is the dipole Hamiltonian that is second order in δM, and H d (r) is the demagnetization field defined by For the infinitely long waveguide, we have H d (r) = 0. For the finite length magnetic bar structure, we approximate H d (r) ≈ H z d (r)ẑ as the z-component is dominant compared to the x and y components. Therefore, we obtain Up to the quadratic order in m(r), we obtain Using M s (r) = M s F (r) and writing m(r) = M(r)F (r), we obtain Similarly, the exchange Hamiltonian can be written, using M(r) = M 0 (r) + δM(r), as Up to the quadratic order in m(r), the above equation becomes Using M s (r) = M s F (r), α ex (r) = α ex F (r), and writing m(r) = M(r)F (r), we obtain where the double-dot product is ∇ M : Note that the term ∂ µ F 3 (r) in the second equation gives a delta-functional contribution peaked at the ferromagnet's boundary. Using the totally-free surface spin condition, ∂ µ M = 0 on the ferromagnet's boundary, we obtain (A29) Combining equations (A23), (A24), (A25), and (A29), we obtain where we dropped the constant shift in energy.
Appendix B: Infinitely long ferromagnetic waveguide

Diagonalization of the magnon Hamiltonian
To obtain the magnon dynamics and the magnon spatial profiles for the infinitely long ferromagnetic waveguide (l → ∞), we diagonalize the magnon Hamiltonian Eq. (A30) by expanding M(r) as where κ X n = nπ/d, κ Y m = mπ/w, n, m = 0, 1, · · · , e ± = x ± iŷ, and a k,(n,m) is the complex canonical variable in the new basis. Note that we have m(r) = M(r)F (r) and in the current geometry F (r) = F X (x)F Y (y), where F X (x) = Θ(x)Θ(d−x), F Y (y) = Θ(y)Θ(w−y), and Θ is the Heaviside step function. Recalling M s (r) = M s F (r) and using Eqs. (A8) and (A9), the above expansion cor-responds to the following: which are presented in the main text. After simplification, the magnon Hamiltonian Eq. (A30) becomes, with where Here, H XX k,(n1m1)(n2m2) , H XY k,(n1m1)(n2m2) , H Y X k,(n1m1)(n2m2) , and H Y Y k,(n1m1)(n2m2) are given by where K α is the modified Bessel function of the second kind, ρ = xx + yŷ, and ϕ XY As we consider the case where the thickness d and the width w are small such that the exchange energy difference D ex K 2 k,(n1m1) − D ex K 2 k,(n2m2) [with (n 1 , m 1 ) = (n 2 , m 2 )] is large as compared to the off-diagonal components [elements of A k,(n1m1)(n2m2) or B k,(n1,m1)(n2m2) with (n 1 m 1 ) = (n 2 , m 2 )] of the Hamiltonian, we apply the block-diagonal approximation [50]. Note that we can go beyond the block-diagonal approximation with the procedure using a paraunitary matrix presented in Ref. [54]. Under this block-diagonal approximation, we obtain and we obtain Now we limit our discussion to the subspace with (n, m) = (0, 0) that gives the lowest energy magnon band, for which magnetization dynamics is uniform across x-y plane in the ferromagnet. After promoting the classical complex canonical variables to the quantum creation and annihilation operators via β k,(0,0) → √ β k,(0,0) and β * k,(0,0) → √ β † k,(0,0) , we obtain which is presented in the main text. Here, ω k,(00) is the magnon energy and β k,(00) is the normal mode magnon annihilation operator satisfying [β k,(00) , β † k ′ ,(00) ] = 2πδ(k − k ′ ). For calculating the dispersion relation in the main text, we numerically evaluate Eqs. (B14)-(B17). In the subspace with (n, m) = (0, 0), ψ X n (x) and ψ Y m (y) are constant functions, so the derivatives only act on F X (x) and F Y (y), resulting in the surface integrals and the evaluation is simpler. Beyond the diagonal approximation (n 1 , m 1 ) = (n 2 , m 2 ) made for Eq. (B8), we can diagonalize the full Hamiltonian via the Bogoliubov transformation with the paraunitary matrix [54] after a truncation of large wavenumber modes, which is used in the magnetic bar calculations in Sec. II and Appendix C.

NV-magnon coupling
The coupling strength between magnons and NV centers is obtained by applying the same Bogoliubov trans-formation in the interaction Hamiltonian Eq. (3). Up to the quadratic order in m(r), we obtain In the infinitely long waveguide case, we have H d (r) = 0. Up to the lowest order (linear order) in m(r), we obtain As the NV axis is setn NV =ẑ, the rotating-wave term comes from the perpendicular contribution h ⊥ (r) = h x (r)x + h y (r)ŷ. Using the Bogoliubov transformation (B22), we obtain where ω d = µ 0 ( γ) 2 /( d 3 ) and Here, Γ XX k,nm , Γ XY k,nm , Γ Y X k,nm , and Γ Y Y k,nm are functions of ρ, and they are given by whereφ XY nm = √ dwϕ XY nm is a dimensionless function and ρ − ρ ′ = (ρ − ρ ′ )/|ρ − ρ ′ |. We consider the external field range γH ext < D NV , where NV center's ground state is |g = |S z = 0 and the first excited state is |e = |S z = − . In the NV center's subspace spanned by {|g , |e }, we can write where ω NV = D NV − γH ext , σ z NV = |e e| − |g g|, and we drop a constant shift in energy. We also have S + NV = √ 2 σ − NV and S − NV = √ 2 σ + NV , where σ + NV = |e g|, and σ − NV = |g e|. Under the rotating wave approximation, we obtain Limiting our discussion to the subspace with (n, m) = (0, 0), we obtain which is presented in the main text. Here, g(ρ i , k) is the dimensionless coupling. To calculate the spatial distribution of the dimensionless coupling, we evaluate Eqs. (B34)-(B37) numerically.

Effective NV-NV Hamiltonian
The NV-NV interaction mediated by magnons can be calculated via the Schrieffer-Wolff transformation [57], is the interaction Hamiltonian in the interaction picture, we obtain the following effective Hamiltonian which is related to the linear response theory. This effective Hamiltonian includes the Lamb shift, the Stark shift, and the NV-NV interaction. The NV-NV interaction contribution is, assuming ρ 1 = ρ 2 and writing g(k) = g(ρ i , k), which is presented in the main text. Here, g eff is the effective NV-NV coupling strength. The entangling gate rate presented in Fig. 2(e) is based on the inverse of the time required for the √ iSWAP gate, τ √ iSWAP = π/(4|g eff |), under the interaction Hamiltonian (B44).
To evaluate how good the perturbation is, we consider one NV case and recall the wave function modification in the first order perturbation where |n (0) and E n are the unperturbed eigenstate and eigenenergy. The fraction of the finite magnon-number state contribution in the original ground state |n (0) = |g |0 m is, where |0 m is the magnon vacuum, Under the geometry presented in the red curve in Fig. 2(e), we obtain n (1) 2 ≈ 10 −3 ≪ 1, which indicates the perturbation theory is valid.
To estimate the corresponding cooperativity of the red solid curve in Fig. 2(e), we assume the waveguide has a length l as in [36]. By discretizing the integral dk using the periodic boundary condition and rescaling the creation/annihilation operators viaβ k,(0,0) = β k,(0,0) / √ l to have a correct commutation relation for the discretized modes, [β k,(0,0) ,β † k ′ ,(0,0) ] = δ k,k ′ , the interaction Hamiltonian becomes As we are mostly using magnons with |k| ≈ k min in the virtual-magnon mediated NV-NV coupling, it is reasonable to calculate the equivalent cooperativity with g =ḡ(k min ): Under the geometry presented in the red curve in Fig. 2(e), and using the NV center's coherence time [7] T * 2 = 1 ms and the Gilbert damping parameter of YIG [67] α = 10 −5 , we obtainḡ ≈ 130 kHz and C eq ≈ 3700.

Temperature independence of the effective NV-NV coupling mediated by virtual magnons
Here we show that up to second order in perturbation theory, the NV-NV coupling mediated by the virtual magnons is insensitive to the temperature. For simplicity, here we only consider the case where two NV centers are coupled to a common single k-magnon mode with coupling strength g k for both NV centers, i.e., , and [a k , a † k ] = 1, although the discussion can be generalized to a multimode or a waveguide case. To demonstrate that, we calculate through the transition matrix formalism the rate T |e1g2n k →|g1e2n k from an initial pure state |e 1 g 2 n k (|n k = (a † k ) n k |0 / √ n k ! with n k = 0, 1, 2, · · · ) to the final state |g 1 e 2 n k , where |i represent the whole set of intermediates manybody states and E |i is the energy of the state |i without interaction. The transition is only non-null for |i = |NV states ⊗ |n k ± 1 , yielding T |e1g2n k →|g1e2n k = 1 g 1 e 2 n k | H int |g 1 g 2 n k + 1 g 1 g 2 n k + 1| H int |e 1 g 2 n k E |g1e2n k − E |g1g2n k +1 (B54) thus first proving the insensitivity to the initial magnon state |n k .
Moreover, we recall that for finite temperature we do not have the pure initial state |e 1 g 2 n k for a specific n k but rather a statistic mix of them, given by the quantum thermal state ρ 0 = Z −1 n k e −β n k ω k |e 1 g 2 n k e 1 g 2 n k |, Z = n k e −β n k ω k with the inverse temperature β = 1/k B T and ω k = ω NV + ∆ k . Finally, using the linearity of the quantum evolution it is straightforward to prove the temperature independence of the off-resonance transition |e 1 g 2 → |g 1 e 2 .
Appendix C: Finite length ferromagnetic bar
1. Firstly, we decomposeĤ into a product of an upper triangle matrix K and its Hermitian conjugate using the Cholesky decomposition 2. Next, we define a new Hermitian matrix W = Kσ 3 K † and diagonalize this matrix with a unitary matrix U: Note that one can find U such that the right-hand side becomes the desired form, which is proven in Ref. [54].
3. Lastly, we define the following matrix T: Then the desired paraunitary matrix is To obtain the eigenfrequencies of the magnons for the finite magnetic bar case, we restrict our discussion for (n, m) = (0, 0) and consider p = 0, 1, · · · , N , where p = N is the highest z-directional wavenumber to be taken into account and we truncated the sum. We set µ 0 = (000), µ 1 = (001), · · · , µ N = (00N ). After the above Bogoliubov transformation with the paraunitary matrix, we obtain with corresponding transformation given by  . . .

NV-magnon coupling
The coupling strength between magnons and NV centers is obtained by applying the same Bogoliubov transformation with the paraunitary matrix T [Eq. (B26)]. Although the demagnetization field H d contribution in (B26) is not negligible when NV centers are placed near the two edges of the ferromagnetic bar, we verify it is small in the calculations for Figs. 4(d) and 5. In the same way as in Sec. I, the perpendicular component of the fringing field h ⊥ is given by where ω dwl = µ 0 (γ ) 2 /( wld) and Here Γ XX µ , Γ XY µ , Γ Y X µ , and Γ Y Y µ are functions of r, and they are given by is a dimensionless function. In the same way as in Sec. I, and under the rotatingwave approximation, we obtain the NV-magnon interaction Hamiltonian in the form of the Jaynes-Cummings which is presented in the main text. To calculate the spatial distribution of the dimensionless coupling, we evaluate numerically Eqs. (C35)-(C38).

Effective NV-NV Hamiltonian
When we introduce a detuning between the target mode frequency ω (00p) and the NV frequency ω NV , we obtain an effective Hamiltonian in the same way as in Sec. I. Now the total Hamiltonian H = H 0 + H int with H 0 = H NV + H m is given by Eqs. (B38), (C29), and (C39). For the the Schrieffer-Wolff transformation, we in the same way as in Eq. (B42). Following Eq. (B43), we obtain where the first right hand side term is the Lamb shift and the Stark shift, respectively. The interaction Hamiltonian between the two NV centers is given by the last right hand side term. If we detune the NV frequency from the mode frequency for µ = (00p) by ω NV = ω (00p) − ∆ω and if we only consider the effect from the mode µ = (00p), we obtain In Fig. 4(c), we focus on the magnon mode with p = 5 and plot the bare coupling g (00p) (r). In Fig. 4(d), we use Eq. (C44) focusing on the magnon mode with p = 5 and plot the effective coupling strength g eff .
Appendix D: Transduction and virtual-magnon exchange protocols

Governing equations for numerical simulations
The comparison between the two entanglement protocols discussed in the main manuscript is performed with the Lindblad master equation simulation [71,72] focusing only on the magnon mode with µ = (00p), p = 5, as presented in Eq. (17). The total Hamiltonian H = H NV + H m + H int to be used is given by Eqs. (B38), (C29), and (C39), where we only considered a single magnon mode µ = (00p) = (005). The identification κ = αω µ presented in Sec. V is appropriate as the dissipation term in the LLG equation ∂ t M| diss = +(α/M s )M × ∂ t M results in ∂ t β µ ≈ −iω µ β µ − αω µ β µ , which is consistent with the master equation result ∂ t a = −iω µ a − κ a when considering only the Boson Hamiltonian. More specifically, as we are considering the case where the equilibrium magnetization is along the z-axis, the linearized equation of motion yields ∂ t m(r) = ∂ t m(r)| coh + αẑ × ∂ t m(r), where ∂ t m(r)| coh is the coherent evolution part described by Eq. (A6). This leads to The positive frequency solutions [solutions with ∼ exp(−iωt)] are obtained by finding nontrivial solutions of Here we notice that one can write H − iαω dra * (r)a(r) = H| ωH →ωH −iαω [See Eqs. (B8)-(B9) and Eqs. (C5)-(C8)], i.e., the Gilbert damping term can be included in the external magnetic field contribution via ω H → ω H − iαω [51,91]. To find ω that gives nontrivial solution, we firstly set α = 0 and obtain ω = ω µ . Then we obtain the solution in the case α = 0 (α ≪ 1) as [51,91] As we can see from Fig. 4 In the simulation presented in Fig. 5, the two NV centers are placed at (x 1 , y 1 , z 1 ) = (d+ h, w, 400 nm) and (x 2 , y 2 , z 2 ) = (d + h, w, 400 nm + δz) with δz = 2.2 µm, which results in g (005) (r 1 ) = g and g (005) (r 2 ) = −g with g = 2π × 517 kHz. The simulation is performed under the field H c , which gives the magnon frequency ω (005) ≈ 2π×2.78 GHz. Moreover, we solve the Lindblad equation in the rotating frame with frequency ω (005) for the transduction protocol and with frequency ω NV for the virtual-magnon exchange protocol. As the NV center's longitudinal relaxation time T 1 is longer than both T * 2 and 1/(αω m ), we do not include its corresponding terms D[σ − NVi ] and D[σ + NVi ] in the current simulation. As shown in the left schematic of Fig. 5(a), idler frequencies of NV 1 and NV 2 in the transduction protocol are ω NV 1 = ω m + δω idle and ω NV 2 = ω m − δω idle , respectively. The detuning δω idle = 2π × 5 MHz is chosen as the neighboring frequencies around ω (005) are separated by more than 2π × 10 MHz from ω (005) , as shown in the Fig. 4(b). The iSWAP gate time is τ iSWAP = π/(2g). Starting from the initial state |g 1 |e 2 , the fidelity is calculated as the state overlap between the NV state and the expected entangled state |ψ ∝ 1 √ 2 (|g 1 |e 2 + e −iδω idle τiSWAP |e 1 |g 2 ). On the other hand, the detuning in the virtual-magnon exchange protocol is ω NV = ω m − ∆ω with ∆ω = 2π × 3 MHz, and the fidelity is calculated as the state overlap with |ψ = 1 √ 2 (|g 1 |e 2 − i|e 1 |g 2 ). The indicator of the violation of the Bell inequality presented in Fig. 5 is calculated following Refs. [73] and [74] as where h j (j = 1, 2, 3) are eigenvalues of the ma- . When (CHSH violation) > 0, the Clauser-Horne-Shimony-Holt (CHSH) form of Bell inequality is violated. As shown in Fig. 5, this is stricter condition than the inseparability of the two-qubit state captured by the entanglement negativity [59], N > 0.

Supplementary simulations
In Fig. 6, we show the temperature dependence of the two entanglement protocols as mentioned in the main text. While we only present the case with T = 70 mK case in Fig. 5, here we present simulations under T = 30 mK, 70 mK, 150 mK, and 300 mK. As the virtualmagnon exchange protocol does not populate the magnon level in the limit ∆ω/g → ∞, i.e., magnons are only created virtually, it is observed that this protocol is robust against the thermal fluctuations. At the same time, as shown in the simulation under T = 30 mK, transduction protocols improves drastically from T = 70 mK compared to the virtual-magnon exchange protocol.
To explore the parameters α and T * 2 dependence of the fidelity on the final entangled state for each protocols, we show in Fig. 7 the parameter dependence of the fidelity at T = 0. The rightmost figure in Fig. 7 shows the phase diagram for which protocol gives better fidelity, where maximum fidelity from each protocols are compared. In the virtual-magnon exchange protocol denoted as detuned, we choose ∆ω = 10g. To simplify the numerical calculation, fidelity at times t = (integer) × π √ 2+(∆ω/g) 2 /g are evaluated for the virtual-magnon exchange protocol, which gives approximately optimal fidelity (see small oscillations observed in the ∆f = 3 MHz cases in Fig. 6). For the transduction protocol, fidelity is evaluated at the time after τ iSWAP /2 interaction time of entangling NV 2 and magnons followed by τ iSWAP iSWAP-gate time between NV 1 and magnons. Here, the coupling strength g µ (r i ) is controlled to be g µ (r i ) = 0 for non-interacting duration instead of inserting idling frequency δω idle , for simplicity. As the resulting fidelity in the virtual-magnon exchange protocol depends on the amount of the detuning ∆ω/g, we show in Fig. 8 the same simulation as in Fig. 7 under multiple detuning values. As shown in the right-top figure in Fig. 8, when the detuning is large ∆ω/g = 30, higher fidelity entangled state can be created even when the magnon damping αω is not very small. This is because magnons are only excited virtually in the virtual-magnon exchange protocol.
As indicated from the phase diagrams presented in Figs. 7 and 8, in the regions where α and T * −1 2 are both sufficiently small, the transduction protocol is better when αω is much smaller than T * −1

2
. On the other hand, virtual-magnon exchange protocol is better when T * −1 2 is much smaller than αω. This tradeoff comes from the fact that the transduction protocol is the faster in gate operation but populate real magnons that are sensitive to the magnon damping, while virtual-magnon exchange protocol is slower in gate operation but it does not populate magnon states and hence the protocol is insensitive to the magnon damping. In Fig. 9, we present the behavior of the boundary line between the two regions for small α and T * −1 2 , where the boundary can be approximated to αω/g = (slope) × (T * −1 2 /g) + (offset). We note that the offset has nodes for detuning values This comes from the small and fast oscillation on top of the slow envelope oscillation observed in the virtual magnon exchange protocol of Fig. 6. The virtual-magnon exchange protocol without the magnon damping and the NV decoherence gives a perfect entangled state only when the condition represented by Eq. (D11) is satisfied. Under this condition, fidelity in the region αω/g ≪ 1 and On the other hand, fidelity in the transduction protocol in the region αω/g ≪ 1 and T * −1 2 /g ≪ 1 is calculated as . Contours indicate Fidelity = 0.5, 0.6, 0.7, 0.8, and 0.9. A phase diagram for which protocol gives better fidelity is presented on the rightmost figure, where the red cross marker represents the parameters used in Fig. 5. We choose ∆ω = 10g for this simulation. For the simplicity of the numerical simulation, we turn on and off the coupling strength instead of inserting the idling frequency δω idle .
Combining Eqs. (D12) and (D13), we obtain the slope value of the boundary line shown in Fig. 9 for detunings ∆ω/g that give zero offset. When ∆ω/g is large, the asymptotic behavior of the slope is (slope) ∼ π 4(π − 1) (∆ω/g) ≈ 0.367(∆ω/g), (D14) which matches with the numerical simulation presented in Fig. 9. However, note that in real magnonic system the detuning is limited by the neighboring mode's frequency separation.
Based on the simulation in Fig. 9, the boundary line under the detuning ∆ω = 2π × 3 MHz is numerically obtained as (αω/g) = 1.24 × 10 −4 + 1.95(T * −1 2 /g). The Gilbert damping parameter α that makes the two protocol comparable is α = 1.35×10 −7 . In Fig. 10, we show the same simulation as in Fig. 5 with parameters T = 0 and α = 1.35 × 10 −7 where we see comparable entanglement values for both protocols, although the transduction protocol is faster in gate operation. For consistency with the analysis presented in Figs. 7-9, the coupling strength g was turned on and off as a function of time instead of inserting the idling frequency δω idle .
Appendix E: Magnon-originated NV center decoherence

Higher order magnon contribution
In this section, we will estimate the decay and decoherence of NV centers due to the interaction with magnon modes with µ = (005) at field H c , which were not taken into account in the Lindblad simulation in the main text. Based on the interaction Hamiltonian Eq. (C39), as the modes with µ = (005) are well separated in frequency, they do not affect the decay and decoherence of NV centers as long as the linewidth αω µ is small. Here we go beyond the linear order interaction, and consider the following NV-magnon interaction (see Eq. (B26)), where h(r) is provided in Eq. (B28) and we define Here, the average is taken with the magnon thermal state ρ m = exp[− µ ω µ β † µ β µ /k B T ], i.e., · · · = Tr[· · · ρ m ]. In the NV center's subspace spanned by {|g , |e }, we can write Assuming a Markovian magnon bath, the NV center's longitudinal decay rates (1/T 1 ) are Under the same assumption, the NV center's decoherence rate (1/T * 2 ) is related to the ω ≈ 0 region of S(ω) with where the Ramsey decoherence follows The longitudinal relaxation rate Γ 1 |e →|g The former is the onemagnon decay contribution (ω NV = ω µ ) and the latter is the two-magnon decay contribution (ω NV = ω µ − ω ν ). However, in our discretized magnon modes, the chances of having ω NV = ω µ or ω NV = ω µ − ω ν are small, at least when the linewidth αω µ of magnons is narrow.
In contrast, for the decoherence that is obtained from ω ≈ 0 part of S(ω), there is a big contribution from terms of the form dte −iωt δn µ (t)δn µ (0) , where δn µ = β † µ β µ − β † µ β µ . This arises from the second-order noise correlation of δh z 2 (r). Furthermore, we notice that this noise contribution is coming not only from the magnon mode with ω µ ≈ ω NV , but also from high energy magnons up to ω µ < k B T / . As the decoherence contribution is expected to be dominant, we estimate the order of its timescale. To simplify the calculation and to avoid the paraunitary matrix diagonalization of a large matrix, we approximate that a µ is the normal mode, i.e. a µ ∼ e −iωµt . We take ω µ = ω min + DK 2 µ , where ω min is the minimum frequency obtained from the paraunitary matrix diagonalization in Sec. II. Hence we write The terms that affect the NV center's decoherence are the contributions from µ = µ ′ . Thus, to estimate the decoherence rate, we take In the limit α → 0 (although this is not compatible with the Markov approximation) we have δn µ (t)δn µ = n 2 µ − n µ 2 that yields where τ 2 is the decoherence timescale. This expression is acceptable as long as the magnon damping 2αω µ is much smaller than 1/τ 2 . When 2αω µ is not small, we take δn µ (t)δn µ = ( n 2 µ − n µ 2 )e −2αωµt and obtain where T * 2 is the decoherence rate. In Fig. 11, we show the two decoherence times from Eqs. (E17) and (E20).

Dispersive coupling contribution
While the Hamiltonian Eq. (C39) does not appear to cause a decoherence, after performing the Schrieffer-Wolff transformation in the dispersive regime (|ω µ − ω NV | > |g µ (r NV )|), we obtain Eq. (C44), where we can securely affirm that the second term (Stark shift term) will cause the decoherence, as considered in Ref. [34]. In this section we calculate the decoherence due to this contribution. We consider the effect of We exclude µ = (005) in the sum as we are considering the field H c where ω NV and ω (005) are on resonant. In the same way as in Eqs. (E17) and (E20), we obtain In Fig. 12, we show the two decoherence times from Eqs.(E23) and (E24). From T ≤ 70 mK and α = 10 −5 part of Figs. 11 and 12, the magnon induced decoherence time is T * 2 > 20 µs, and it is expected that this dephasing contribution does not change the general trend of the result of the simulation presented in Fig. 5.
In Fig. 13, we show figures corresponding to Fig. 6 with the NV centers' dephasing rate calculated in Fig. 12. The simulation confirms that the general tendency presented in Fig. 5 does not change due to the dephasing contribution calculated in Fig. 12.
Appendix F: Average Gate Fidelity for off-resonance protocol To show that the magnon-mediated entanglement protocols can directly be extended to two-qubit gates, in this section we have calculated the average gate fidelity as a square-root-of-iSWAP gate for the off-resonant protocol under the same condition as in Fig. 5. To calculate the average gate fidelity, we employ a method based on the entanglement fidelity F e [75]. For that we introduce two auxiliary qubits aux 1 and aux 2 and prepare the following maximally entangled state [92] as an initial qubit state. Then we evolve in time the NV and magnon states according to the Lindblad master equation of the previous sections, and calculate the fidelity F e as the state overlap between the calculated state and the desired state after the following gate where τ √ iSWAP = π/(4|g eff |). As the square of U gate is equivalent to the iSWAP gate up to single-qubit operations, U gate can be thought of as a square-root-of-iSWAP gate. The average gate fidelityF is calculated via [75]  The Gilbert damping parameter α = 10 −5 that is observed in bulk YIG crystals [67] would be optimistic for small YIG structures that we consider in this work. However, as one can calculate from Fig. 4(c), we obtain a high cooperativity C ≈ 500 even with a larger Gilbert damping parameter α = 10 −3 . In Fig. 15, we show a simulation analogous to the one presented in Fig. 5 with α = 10 −3 . From this simulation, we find that the off-resonance protocol produces entangled states, as the entanglement negativity is larger than zero. However, this turns out to be not a useful entanglement as (CHSH Violation) = 0 indicates that the state does not violate the Bell inequality. This happens because of the increased T 1 decay rate of NV centers due to the overlap of the broad magnon mode resonance with the NV-center's transition. Although the off-resonance protocol is less sensitive to the magnon decay, the detuning ∆ω needs to be sufficiently larger than the linewidth of the magnon-mode resonance αω µ in order to suppress this decay contribution.
The resulting entangled mixed state presented in Fig. S10 can be understood in the following way. As the interaction Hamiltonian is H int = ga(σ + NV1 −σ + NV2 )+H.c., we notice that |D = (|g NV1 |e NV2 + |e NV1 |g NV2 )/ √ 2 is a dark state with respect to the magnon mode, or alternatively, |D is a state within a subspace that is free from the magnon-induced T 1 decay (decoherence free subspace), because H int |D |n m = 0 with a magnon number state |n m . Accordingly, the initial state of NV cen- ters can be written as |ψ init = |g NV 1 |e NV 2 = (|D + |B )/ √ 2, with |B = (|g NV1 |e NV2 − |e NV1 |g NV2 )/ √ 2, and initial density operator ρ init = |ψ init ψ init | = (|D D| + |D B| + |B D| + |B B|)/2. After the time evolution, the part related to |D D| remains constant as |D is in the decoherence free subspace. Assuming that the system is at absolute zero temperature for simplicity, and that the other terms eventually evolve to the ground state |D B| + |B D| + |B B| → |00 00| due to the energy relaxation, where |00 = |g NV1 |g NV2 , we obtain the final density operator ρ fin = (|D D| + |00 00|)/2. (G1) As the partial transpose of this density matrix has a negative eigenvalue −( √ 2 − 1)/4, we obtain the entanglement negativity of the final state N fin = ( √ 2 − 1)/4 and N fin /N B = ( √ 2 − 1)/2 ≈ 0.21. This explains the lower-right panel of Fig. 15 with an additional note that at T = 70 mK the final density operator that evolved from |D B| + |B D| + |B B| is no longer |00 00|, but rather a mixture of |00 00|, |B B|, and |11 11|, where |11 = |e NV 1 |e NV 2 .
To mitigate the magnon-induced T 1 decay in the case of the larger Gilbert damping parameter, one can make the detuning ∆f larger. Although in our case this is limited by the frequency spacing of the neighboring magnon modes [see Fig. 4(b)], we show in Fig. 16 the simulation with a larger detuning value ∆f = 30 MHz. We note, however, that this is not possible for the magnonic system we have considered in the main text as the neighboring magnon-mode frequency separations are smaller than 30 MHz [See Fig. 4(b)] in the main text. Conversely, this simulation clarifies that the system will make useful entanglement that can violate the Bell inequality. This implies that to improve the quality of the resulting entanglement further optimization on the length l of the magnetic bar structure is needed, as it defines the frequency spacing of magnon modes. β † µ (0)β v (t) = β † µ (0)β v (0) e −iωµt−|κ|t with κ = αω µ , we obtain [93,94] γ 1 −,X = µ=(00p) p=0,1,··· g (X) p 2 (n B (ω µ ) + 1) · 2κ (Ω X − ω µ ) 2 + κ 2 , ≈ µ=(00p) p=0,1,··· g (X) p 2 (n B (ω µ ) + 1) · 2κ (Ω X − ω µ ) (Ω X − ω µ ) 2 + κ 2 , ≈ µ=(00p) p=0,1,··· g (X) p 2 n B (ω µ ) · 2κ (Ω X − ω µ ) 2 , where n B (ω) = [exp( ω/k B T )− 1] −1 is the Bose-Einstein distribution function and we have approximated (Ω L/U − ω µ ) ≫ κ to obtain the last expressions. Note that for the lower frequency transition, we do not include p = 5 in the summation as this is the on-resonant magnon mode and its effect is directly included in the simulation in Fig. 5. With the Gilbert damping parameter α = 10 −5 , we evaluated the above expression and obtained Fig. 18. As the calculated relaxation time is much longer than the time scale that is simulated in Fig. 5, this T 1 decay contribution from magnon modes other than p = 5 is negligible for the condition we considered.
According to Eqs. (B43) and (B45), the effective NV-NV interaction is related to the ω = ω NV contribution of the Green's function G A (ω NV ). Assuming D k,k ′ will contribute to the effective NV-NV coupling on the same order as 2πδ(k − k ′ ) in Eq. (K12) for simplicity to evaluate the scale of the contribution of the perturbation and as they have the same dimension of length, and using ω h G 0 (ω NV , k ′ ) ∼ ω h /(ω NV − ω k ′ ,(0,0) ) ∼ ω h /(ω min − ω NV ), the effect of the local magnetic field h ext on the NV-NV effective coupling, based on Eqs. (B43),(B45), and (K12), is given by although further investigation is needed for the full comparison of the two terms in Eq. (K12) as well as for higher order terms.

Perturbative approach to the YIG bar case
In the case of the YIG bar, with the use of Eqs. (C1),(C2), the Hamiltonian H (2) m can be written as Now we define the perturbation Hamiltonian matrix λV as In the following, we will consider the effect of λV in the expansion with the order λ for the case of the diagonalization with a paraunitary matrix. We want to diagonalize the total Hamiltonian matrixĤ =Ĥ 0 + λV in the form and we assume we know this expansion in the case with λ = 0 as Based on these, we expand the perturbed paraunitary T and eigenvalues Λ matrices as T = T 0 + λT 1 + · · · . (K21) Λ = Λ 0 + λΛ 1 + · · · .
(K22) Substituting these into Eqs. (K17) and (K18), and taking leading order terms in λ, we obtain where D is an arbitrary diagonal matrix with purelyimaginary entries. This is due to the degrees of freedom of the paraunitary matrix T → Texp[iλ × (real diagonal matrix)], which we encounter in the unitary diagonalization case as well. Therefore, we simply set D = 0 and obtain where g 0 (00p) is the coupling strength we obtained without the perturbation Hamiltonian H (2) m . From Eqs. (K23) and (K26) with Λ 1 ∼ ω h and L ∼ ω h /(ω ν − ω µ ), we find the following scaling behavior for the change in the magnon mode frequency and the NV-magnon coupling strength due to the local nonuniform magnetic field h ext , although Eq. (K32) strongly depends on how much the additional magnetic field mixes different normal magnon modes, described by the off-diagonal components of T † 0 VT 0 .