Quantum entanglement between two magnon modes via Kerr nonlinearity

TL;DR

This paper proposes a method to generate and maintain quantum entanglement between two macroscopic magnon modes in YIGs using Kerr nonlinearity and strong microwave driving, with potential applications in quantum technologies.

Contribution

It introduces a scheme to entangle two ferromagnetic samples via Kerr nonlinearity in a driven cavity system, demonstrating steady-state entanglement in massive magnetic materials.

Findings

Entanglement between two magnon modes is achieved via Kerr nonlinearity.

Steady-state entanglement persists under strong microwave driving.

The scheme provides a pathway for observing quantum effects in macroscopic magnetic systems.

Abstract

We propose a scheme to entangle two magnon modes via Kerr nonlinear effect when driving the systems far-from-equilibrium. We consider two macroscopic yttrium iron garnets (YIGs) interacting with a single-mode microcavity through the magnetic dipole coupling. The Kittel mode describing the collective excitations of large number of spins are excited through driving cavity with a strong microwave field. We demonstrate how the Kerr nonlineraity creates the entangled quantum states between the two macroscopic ferromagnetic samples, when the microcavity is strongly driven by a blue-detuned microwave field. Such quantum entanglement survives at the steady state. Our work offers new insights and guidance to designate the experiments for observing the entanglement in massive ferromagnetic materials. It can also find broad applications in macroscopic quantum effects and magnetic spintronics.

Equations

$$\begin{array}{rl}H=& -\int M_z{B}_0\text{d}\text{r}-\frac{{\mu }_0}{2}\int M_z{H}_{\text{an}}\text{d}\text{r}\\ & +\frac{1}{2}\int \left({\epsilon }_0{\text{E}}^2+\frac{{\text{B}}^2}{{\mu }_0}\right)\text{d}\text{r}-\int \text{M}\cdot \text{B}\text{d}\text{r}\end{array}$$

$$\begin{array}{rl}H=-\gamma & \sum _{j=1}^2{B}_{j,0}{S}_{j,z}+{\gamma }^2\sum _{j=1}^2\frac{{\mu }_0{K}_{\text{an}}^{(j)}}{M_j^2{V}_{j,m}}{S}_{j,z}^2\\ & +\hslash {\omega }_c{a}^{†}a-\gamma \sum _{j=1}^2{S}_{j,x}{B}_{j,x}\end{array}$$

$$\begin{split}H_{\text{eff}}=\hbar&\omega_{c}a^{\dagger}a+\hbar\sum_{j=1}^{2}\Big{[}\omega_{j}m_{j}^{\dagger}m_{j}+g_{j}(m_{j}^{\dagger}a+m_{j}a^{\dagger})\\[4.26773pt] &+\Delta_{j}m_{j}^{\dagger}m_{j}m_{j}^{\dagger}m_{j}\Big{]}+i\hbar\Omega(a^{\dagger}e^{-i\omega_{\text{d}}t}-ae^{i\omega_{\text{d}}t})\end{split}$$

$$\begin{array}{rl}& {\dot{m}}_s=-(i{\delta }_s+{\gamma }_s)m_s-2i{\Delta }_s{m}_s^{†}{m}_s{m}_s-i{g}_{s}a+\sqrt{2{\gamma }_s}{m}_s^{\text{in}}(t)\\ & \dot{a}=-(i{\delta }_c+{\gamma }_c)a-i\sum _{j=1}^2{g}_j{m}_j+\Omega +\sqrt{2{\gamma }_c}{a}^{\text{in}}(t)\end{array}$$

$$\begin{array}{r}\dot{\sigma }(t)=A\sigma (t)+f(t)\end{array}$$

$$\begin{array}{r}A=\left(\begin{array}{cccccc}{F}_1-{\gamma }_1& {\tilde{\delta }}_1-G_1& 0& 0& 0& g_1\\ -{\tilde{\delta }}_1-G_1& -F_1-{\gamma }_1& 0& 0& -g_1& 0\\ 0& 0& F_2-{\gamma }_2& {\tilde{\delta }}_2-G_2& 0& g_2\\ 0& 0& -{\tilde{\delta }}_2-G_2& -F_2-{\gamma }_2& -g_2& 0\\ 0& g_1& 0& g_2& -{\gamma }_c& {\delta }_c\\ -g_1& 0& -g_2& 0& -{\delta }_c& -{\gamma }_c\end{array}\right)\end{array}$$

$$\begin{array}{rl}& ⟨m_1⟩=\frac{i{g}_1\Omega }{({\tilde{\delta }}_1-i{\gamma }_1)({\delta }_c-i{\gamma }_c)-g_1^2-\frac{g_2^2({\tilde{\delta }}_1-i{\gamma }_1)}{{\tilde{\delta }}_2-i{\gamma }_2}},\\ & \text{and}(1\leftrightarrow 2)\end{array}$$

$$\begin{split}H_{\text{qf}}=\hbar\sum_{s=1}^{2}&\Big{[}\tilde{\delta}_{s}\delta m_{s}^{\dagger}\delta m_{s}+\tilde{\Delta}_{s}\delta m_{s}^{\dagger}\delta m_{s}^{\dagger}+\tilde{\Delta}_{s}^{*}\delta m_{s}\delta m_{s}\\[4.26773pt] &+g_{s}\left(\delta m_{s}^{\dagger}\delta a+\delta m_{s}\delta a^{\dagger}\right)\Big{]}+\hbar\delta_{c}\delta a^{\dagger}\delta a\end{split}$$

$$\begin{split}H_{\text{qf}}=\hbar&\sum_{s=1}^{2}\Big{[}\varepsilon_{s}\delta\beta_{s}^{\dagger}\delta\beta_{s}+g_{s}\Big{(}(v_{s}\delta\beta_{s}+u_{s}\delta\beta_{s}^{\dagger})\delta a\\[4.26773pt] &\quad+(u_{s}^{*}\delta\beta_{s}+v_{s}^{*}\delta\beta_{s}^{\dagger})\delta a^{\dagger}\Big{)}\Big{]}+\hbar\delta_{c}\delta a^{\dagger}\delta a\end{split}$$

$$\begin{array}{r}\dot{C}(t+\tau ,t)=AC(t+\tau ,t)+C(t+\tau ,t)A^{\text{T}}+e^{A\tau }D\end{array}$$

$$\begin{array}{r}A{C}_{\infty }+C_{\infty }{A}^{\text{T}}=-D\end{array}$$

Full text

Introduction.– Recent advance in ferromagnetic materials draw considerable attention in the studies of quantum nature in magnetic systems, as the limitations of electrical circuitry are reached. Thanks to the low loss of the collective excitations of spins known as magnons in magnetic samples, the magnons offer a new paradigm for developing future generation of spintronic devices and quantum engineering Kajiwara_Nature2010 ; Cornelissen_NatPhys2015 ; Zhu_APL2016 ; An_NatMater2013 ; Chumak_NatPhys2015 ; Chumak_NatCommun2014 . The yttrium iron garnet (YIG) with the size of $\sim 100\mu$m as fabricated in recent experiments provides new insights for studying the macroscopic quantum effects, such as entanglement and squeezing that have raised widespread interest in different branches of physics during decade Collet_NatCommun2016 ; Ho_PRL2018 ; Simon_Nat2018 ; Yuan_PRB2018 ; Morimae_PRA2005 ; Korppi_Nature2018 . Quantum entanglement between massive mirror and optical cavity photons has been explored, in both theoretical and experimental aspects Vitali_PRL2007 ; Simon_Nature2009 ; Aspelmeyer_RMP2014 ; Genes_PRA2008 ; Verhagen_Nature2012 ; Palomaki_Sci2013 . Several ideas follow-on suggest the extension of such entangled quantum state towards the magnons in microwave regime, due to their great potential for macroscopic spintronic devices. Much experimental efforts have been devoted to the quantum nature of magnon states, through hybridizing the spin waves with other degrees of freedoms, e.g., superconducting qubits and phonon modes Julsgaad_Nat2001 ; Wellstood_Science2003 ; Quirion_SciAdv2017 ; Zhang_SciAdv2016 . Compared to atoms and photonics, magnonics holds the potential for implementing quantum states in more massive objects. This can be seen from the $320\mu$m-diam YIG spheres implemented in recent experiments Zhang_NPJ2015 .

As a powerful platform for investigating the light-matter interaction Zhang_NPJ2015 ; Tabuchi_PRL2014 ; Soykal_PRL2010 ; Zhang_PRL2014 ; Bourhill_PRB2016 ; Tabuchi_Sci2015 ; Harder_PRL2018 ; Yao_NatCommun2017 , ferromagnetic materials are taking the advantage of reaching strong and ultrastrong coupling regimes, along with the fact of their high spin density as well as low dissipation rate. The strong coupling results in the cavity magnon-polariton, serving as a potential candidate for implementing quantum information transducers and memories Yao_NatCommun2017 ; Zhang_NatCommun2015 . To achieve the quantum regime in magnon polaritons, the macroscopic quantum effects are essentially worthy of being explored. The most recent work using driven-dissipation theory suggest the magnon-photon-phonon entanglement and also the squeezing of magnon modes in which both the entanglement and squeezing are essentially transferred into the mechanical mode Li_PRL2018 ; Li_PRA2019 ; Li_magnon2019 . From a theoretical view-point, this macroscopic quantum nature of magnon modes stems from the nonlinearity that can be enhanced by driving the systems far-from-equilibrium. Two prominent schemes are responsible for introducing such nonlinearity: the magnetostrictive interaction and the Kerr effect, where the latter results from the magnetocrystalline anisotropy. Apart from the magnon-phonon interaction, Kerr nonlinearity plays a significant role in magnon spintronics Chumak_NatPhys2015 . Recent experiments in YIG spheres demonstrated the multistability and photon-mediated control of spin current, due to the Kerr effect You_PRL2018 ; Hu_PRL2017 .

In this Letter, we propose a scheme of entangling magnon modes in two massive YIG spheres via the Kerr nonlinearity. The two magnon modes interact with a microcavity through the beam-splitter-like coupling, which cannot produce any entanglement. Nevertheless, activating the Kerr nonlinearity via strong driving results in squeezing-like coupling which may let magnon get entangled with cavity photons. The subsequent entanglement transfer between photons and the other magnon mode will lead to the entanglement between magnon modes. The condition for optimizing the magnon-magnon entanglement is found and is confirmed by our numerical calculations. By taking into account the experimentally feasible parameters, we show the considerable magnon-magnon entanglement can be created. Such entanglement is also shown to be robust against cavity leakage. Our work offers new insight and perspective for studying the quantum effects in complex molecules. These have been manifested by the excited-state dynamics in dye molecules and even bacterias implying the entangled quantum states when interacting with microcavities Sarovar_NatPhys2010 ; Zhang_SR2016 ; Coles_NatMater2014 ; Zhang_JCP2018 ; Mukamel_JPCL2016 .

Model and equation of motion.– We consider a hybrid magnon-cavity system consisting of two bulk ferromagnetic materials and one microwave cavity mode. The ferromagnetic sample contains dispersive spin waves, in which only the spatially uniform mode (Kittel mode Kittel_PR1948 ) is assumed to strongly interact with cavity photons. The full Hamiltonian of this cavity magnonics system reads Blundell_book2001

[TABLE]

where $\textbf{B}_{0}=B_{0}\textbf{e}_{z}$ is the applied static magnetic field and $\textbf{M}=\gamma\textbf{S}/V_{m}$ with $\gamma=e/m_{e}$ denoting the gyromagnetic ratio. S stands for the collective spin operator and $V_{m}$ is volume of ferromagnetic material. $\textbf{H}_{\text{an}}$ is the anisotropic field due to the magnetocrystalline anisotropy and has $z$ component only owing to the crystallographic axis being aligned along the applied static magnetic field. Thereby the anisotropic field is given by $H_{\text{an}}=-2K_{\text{an}}M_{z}/M^{2}$ where $K_{\text{an}}$ and $M$ denote the dominant 1st anisotropy constant and the saturation magnetization, respectively. One can recast the Hamiltonian in Eq.(1) into

[TABLE]

by assuming the magnetic field inside cavity is along $x$-axis. The Holstein-Primakoff transform yields to $S_{i,z}=S_{i}-m_{i}^{\dagger}m_{i},\ S_{i,+}=(2S_{i}-m_{i}^{\dagger}m_{i})^{1/2}m_{i},\ S_{i,-}=m_{i}^{\dagger}(2S_{i}-m_{i}^{\dagger}m_{i})^{1/2}$ where $S_{i,\pm}\equiv S_{i,x}\pm iS_{i,y}$ and $m_{i}$ represents the bosonic annihilation operator Madelung_book1978 . For the yttrium iron garnets (YIGs) with diameter $\text{d}=40\mu$m, the density of ferrum ion $\text{Fe}^{3+}$ is $\rho=4.22\times 10^{27}$m*-3*, which leads to the total spin $S=\frac{5}{2}\rho V_{m}=7.07\times 10^{14}$. This is often much larger than the number of magnons, so that we can safely approximate $S_{j,+}\simeq\sqrt{2S_{j}}m_{j},\ S_{j,-}\simeq\sqrt{2S_{j}}m_{j}^{\dagger}$. In the presence of external microwave driving field, the effective Hamiltonian of hybrid magnon-cavity system is of the form

[TABLE]

where the rotating-wave approximation was employed and cavity frequency is denoted by $\omega_{c}$. The frequency of Kittel mode is $\omega_{j}=\gamma B_{j,0}$ with $\gamma/2\pi=28$GHz/T. $g_{j}$ gives the magnon-cavity coupling and $\Delta_{j}=\mu_{0}K_{\text{an}}^{(j)}\gamma^{2}/M_{j}^{2}V_{j,m}$ gives the Kerr nonlinearity, resulting from the on-site magnon-magnon scattering. The Rabi frequency $\Omega=\sqrt{2P_{\text{d}}\gamma_{c}/\hbar\omega_{\text{d}}}$ in last term quantifies the strength of the field inside microcavity driven by the microwave magnetic field, where $P_{\text{d}}$ and $\omega_{\text{d}}$ represent the power and frequency of the microwave field, respectively. $\gamma_{c}$ is the cavity leaking rate. In the rotating frame of microwave field, the dynamics of hyrid cavity-magnon system is governed by the quantum Langevin equations (QLEs)

[TABLE]

where $\gamma_{s}$ quantifies the magnon dissipation. $\delta_{s}=\omega_{s}-\omega_{\text{d}},\ \delta_{c}=\omega_{c}-\omega_{\text{d}}$. $m_{s}^{\text{in}}(t)$ and $a^{\text{in}}(t)$ are the input noise operators having zero mean and white noise: $\langle m_{s}^{\text{in},\dagger}(t)m_{s}^{\text{in}}(t^{\prime})\rangle=\bar{n}_{s}\delta(t-t^{\prime}),\ \langle m_{s}^{\text{in}}(t)m_{s}^{\text{in},\dagger}(t^{\prime})\rangle=(\bar{n}_{s}+1)\delta(t-t^{\prime});\ \langle a^{\text{in},\dagger}(t)a^{\text{in}}(t^{\prime})\rangle=0,\ \langle a^{\text{in}}(t)a^{\text{in},\dagger}(t^{\prime})\rangle=\delta(t-t^{\prime})$ where $\bar{n}_{s}=[\text{exp}(\hbar\omega_{s}/k_{B}T)-1]^{-1}$ denotes the Planck factor of the $s$-th magnon mode.

Since the microcavity is under strong driving by the microwave field, the beam-splitter-like coupling between magnons and cavity leads to the large amplitudes of both magnon and cavity modes, namely, $|\langle m_{s}\rangle|,\ |\langle a\rangle|\gg 1$. In this case, one can safely introduce the expansion $m_{s}=\langle m_{s}\rangle+\delta m_{s},\ a=\langle a\rangle+\delta a$ in the vicinity of steady state, by neglecting the higher-order fluctuations of the operators. We thereby obtain the linearized QLEs for the quadratures $\delta X_{s},\delta Y_{s},\delta X,\delta Y$ defined as $\delta X_{1}=(\delta m_{1}+\delta m_{1}^{\dagger})/\sqrt{2},\ \delta Y_{1}=(\delta m_{1}-\delta m_{1}^{\dagger})/i\sqrt{2},\ \delta X_{2}=(\delta m_{2}+\delta m_{2}^{\dagger})/\sqrt{2},\ \delta Y_{2}=(\delta m_{2}-\delta m_{2}^{\dagger})/i\sqrt{2},\ \delta X=(\delta a+\delta a^{\dagger})/\sqrt{2},\ \delta Y=(\delta a-\delta a^{\dagger})/i\sqrt{2}$

[TABLE]

where $\sigma(t)=[\delta X_{1}(t),\delta Y_{1}(t),\delta X_{2}(t),\delta Y_{2}(t),\delta X(t),\delta Y(t)]^{\text{T}}$ and $f(t){=}[\sqrt{2\gamma_{1}}X_{1}^{\text{in}}(t),\!\sqrt{2\gamma_{1}}Y_{1}^{\text{in}}(t),\!\sqrt{2\gamma_{2}}X_{2}^{\text{in}}(t),\!\sqrt{2\gamma_{2}}Y_{2}^{\text{in}}(t),\!\sqrt{2\gamma_{c}}X^{\text{in}}(t),\!\\ \sqrt{2\gamma_{1}}Y^{\text{in}}(t)]^{\text{T}}$ are the vectors for quantum fluctuations and noise, respectively. The drift matrix reads

[TABLE]

with magnetocrystalline anisotropy quantified by $G_{s}=2\Delta_{s}\text{Re}\langle m_{s}\rangle^{2},\ F_{s}=2\Delta_{s}\text{Im}\langle m_{s}\rangle^{2}$ and the effective detuning of magnons $\tilde{\delta}_{s}=\delta_{s}+2\sqrt{G_{s}^{2}+F_{s}^{2}}=\delta_{s}+4\Delta_{s}|\langle m_{s}\rangle|^{2}$, which includes the frequency shift caused by Kerr nonlinearity. The mean $\langle m_{1,2}\rangle$ are given by

[TABLE]

Before the study of entanglement, it is essential to elucidate the mechanism for optimizing the entanglement via Kerr nonlinearity. To this end, we proceed via the effective Hamiltonian for quantum fluctuations

[TABLE]

where $\tilde{\Delta}_{s}=(G_{s}+iF_{s})/2$. The quadratic terms $\delta m_{s}^{\dagger}\delta m_{s}^{\dagger},\ \delta m_{s}\delta m_{s}$ imply the effective magnon-magnon interaction induced by the magnetocrystalline anisotropy, which may be significantly enhanced by strong driving. This, in fact, is responsible for the entanglement. To make it elaborate, let us introduce the Bogoliubov transformation Lifshitz_book1980 ; Fetter_book2003 $\delta\beta_{s}=u_{s}\delta m_{s}-v_{s}^{*}\delta m_{s}^{\dagger},\ \delta\beta_{s}^{\dagger}=-v_{s}\delta m_{s}+u_{s}^{*}\delta m_{s}^{\dagger}$ where $u_{s}=\sqrt{\frac{1}{2}\left(\frac{\tilde{\delta_{s}}}{\varepsilon_{s}}+1\right)},\ v_{s}e^{i\alpha}=-\sqrt{\frac{1}{2}\left(\frac{\tilde{\delta_{s}}}{\varepsilon_{s}}-1\right)},\ \alpha=\text{arctan}(F_{s}/G_{s})$ and $\varepsilon_{s}=\left(\tilde{\delta}_{s}^{2}-4|\tilde{\Delta}_{s}|^{2}\right)^{1/2}$. Inserting these into Eq.(8) we find

[TABLE]

which shows $\varepsilon_{s}\simeq-\delta_{c}$ is optimal for the entanglement, due to the magnon-photon squeezing term $g_{s}(v_{s}\delta\beta_{s}\delta a+v_{s}^{*}\delta\beta_{s}^{\dagger}\delta a^{\dagger})$. This will be confirmed by the latter numerical results when taking into account of experimental parameters.

Entanglement between magnon modes.– Since we are using the linearized quantum Langevin equations, the Gaussian nature of the input states will be preserved during the time evolution of systems. The quantum fluctuations are thus the continuous three-mode Gaussian state, which is completely characterized by an $6\times 6$ covariance matrix (CM) defined as $C_{ij}(t,t^{\prime})=\frac{1}{2}\langle\sigma_{i}(t)\sigma_{j}(t^{\prime})+\sigma_{j}(t^{\prime})\sigma_{i}(t)\rangle;\ (i,j=1,2,\cdots,6)$ where the average is taken over the system and bath degrees of freedoms. Suppose the drift matrix $A$ is negatively defined, the solution to Eq.(5) is $\sigma(t)=M(t)\sigma(0)+\int_{0}^{t}M(s)f(t-s)\text{d}s$ where $M(t)=\text{exp}(At)$. This enables us to find the equation which CM obeys

[TABLE]

for $\tau\geq 0$. Thus the stationary CM can be straightforwardly obtained by letting $\tau=0,\ t\rightarrow\infty$ in Eq.(10) that yields to the Lyapunov equation

[TABLE]

where the diffusion matrix is $D=\text{diag}[\gamma_{1}(2\bar{n}_{1}+1),\gamma_{1}(2\bar{n}_{1}+1),\gamma_{2}(2\bar{n}_{2}+1),\gamma_{2}(2\bar{n}_{2}+1),\gamma_{c},\gamma_{c}]$ defined through $\langle f_{i}(t)f_{j}(t^{\prime})+f_{j}(t^{\prime})f_{i}(t)\rangle=2D_{ij}\delta(t-t^{\prime})$. We adopt the logarithmic negativity $E_{N}$ to quantify the magnon-magnon and magnon-photon entanglements by comupting the $4\times 4$ CM related to the two magnon modes. This can be achieved by defining $E_{N}=\text{max}[0,-\text{ln}2v_{-}]$ where $v_{-}=\text{min}|\text{eig}\oplus_{j=1}^{2}(-\sigma_{y})P_{12}C_{\infty}P_{12}|$ and $\sigma_{y}$ is the Pauli matrix Vidal_PRA2002 ; Simon_PRL2000 . The matrix $P_{12}=\sigma_{z}\oplus 1$ realizes the partial transposition at the level of CM. In what follows, we will work in the monostable scheme of magnons. Furthermore, we will focus on the case of two identical magnons having $G_{1,2}=G,\ F_{1,2}=F,\ \tilde{\delta}_{1,2}=\tilde{\delta},\ \Delta_{1,2}=\Delta,\ g_{1,2}=g$.

Fig.2 shows the magnon-magnon entanglement versus some key parameters of the system. Here we have taken into account the experimentally feasible parameters You_PRL2018 : $\omega_{1,2}/2\pi=10$GHz, $\delta_{1,2}/2\pi=-1$MHz, $\gamma_{1,2}/2\pi=8.8$MHz and $\gamma_{c}/2\pi=1.9$MHz for the YIG bulk at low temperature $T=10$mK. First of all we observe from Fig.2(a,b) that the Kerr nonlinearity is responsible for creating the steady-state entanglement between two magnon modes, evident by the fact that the entanglement dies out when $G=F=0$. This results from the dominated beam-splitter-interaction between magnon mode and cavity photons, once $G=F=0$. Thereby no magnon-cavity entanglement can be created, as seen in Fig.2(b). We take the condition $\varepsilon_{s}\simeq-\delta_{c}$ for optimizing the magnon-photon entanglement, as illustrated in Fig.2(d) where $\varepsilon_{s}\simeq\sqrt{3(G_{s}^{2}+F_{s}^{2})}$. The two-mode squeezing term $g_{s}(v_{s}\delta\beta_{s}\delta a+v_{s}^{*}\delta\beta_{s}^{\dagger}\delta a^{\dagger})$ squeezes the joint state between one magnon mode and cavity photons, which results in the partial entanglement in between. Because the same type of interaction occurs when coupling the other magnon mode with cavity, the two distanced magnon modes are expected to be entangled. This is confirmed in Fig.2(c) manifesting the optimal magnon-magnon entanglement in the vicinity of $\varepsilon_{s}\simeq-\delta_{c}$. The elaborate transfer from magnon-photon entanglement to magnon-magnon entanglement is subsequently evident by Fig.1S in supplementary material (SM) that the considerable reduction of magnon-photon entanglement as the coupling of cavity to another sphere is turned on. Since the biparticle entanglement is originated from the Kerr nonlinearity quantified by $G_{1,2}$ and $F_{1,2}$, there must be the interplay between the couplings $G_{s},\ F_{s}$ and $g_{s}$ which is depicted in Fig.2(e,f). In Fig.2(c) we take $g_{1,2}/2\pi=41$MHz and it implies $\delta_{c}/2\pi\simeq-0.03$GHz for the optimal entanglement $E_{m_{1}m_{2}}$. We then adopt the magnitude of $\delta_{c}$ for plotting Fig.2(a,b). Using $\sqrt{G^{2}+F^{2}}=2\Delta|\langle m\rangle|^{2}$ and Eq.(7) for the $40\mu$m-diam YIG spheres, the optimal entanglement with $|G|=0.038\text{GHz},\ |F|=0.028$GHz (see Fig.2(a)) yields to the Rabi frequency $\Omega=1.06\times 10^{15}$Hz, corresponding to the drive power $P_{\text{d}}=314$mW. Indeed, the stronger nonlinearity will create more entanglement between the magnon modes. But we have to ensure the negatively defined matrix $A$ given in Eq.(6). Also, the experimental feasibility of ultrastrong drive using microwave field needs the consideration.

Since we are working with the strong driving, it is worthy of checking the validity of the results obtained above. The magnon description for magnetic materials is effective only when $\langle m_{j}^{\dagger}m_{j}\rangle\ll 2N_{j}s=5N_{j}$ where $N_{j}=\rho_{j}V_{j,m}$ denotes the total number of spins in the bulk material. For the $40\mu$m-diam YIG sphere, $N_{j}\simeq 1.41\times 10^{14}$ and the drive power $P_{\text{d}}=393$mW results in $|\langle m_{j}\rangle|\simeq 2.3\times 10^{6}$, giving $\langle m_{j}^{\dagger}m_{j}\rangle\simeq 5.28\times 10^{12}\ll 5N=7.07\times 10^{14}$. Hence the condition $\langle m_{j}^{\dagger}m_{j}\rangle\ll 2N_{j}s$ is fulfilled.

Fig.3 illustrates the entanglement between two magnon modes versus some controllable parameters by considering the $40\mu$m-diam YIG-sphere experiment, where $\omega_{1,2}/2\pi=10$GHz, $\delta_{1,2}/2\pi=-1$MHz, $\Delta_{1,2}/2\pi=1\mu$Hz, $g_{1,2}/2\pi=41$MHz, $\gamma_{1,2}/2\pi=8.8$MHz and $\gamma_{c}/2\pi=1.9$MHz have been taken according to Ref.You_PRL2018 . We observe in Fig.3(a) that for fixed driving power, the magnon-magnon entanglement is quite sensitive to cavity detuning $\delta_{c}\equiv\omega_{c}-\omega_{\text{d}}$, reaching its maximum at $\delta_{c}/2\pi\simeq-0.03$GHz. This is consistent with the condition $\varepsilon_{j}\simeq-\delta_{c}$ as clarified for optimizing the entanglement. Fig.3(b) shows the considerable entanglement when the system is driven far-from-equilibrium. This is reasonable because the strong external driving significantly enhances the Kerr nonlinearity that is responsible for both magnon-cavity squeezing and entanglement, as elucidated in Eq.(7) and (8). Furthermore, Fig.3 shows that the weaker magnon-magnon entanglement is observed when increasing the cavity leakage. By noting the magnitude, we can still obtain some entanglement, even with a low-quality cavity showing weak magnon-cavity coupling where $\gamma_{c}=8\gamma_{1,2}>g_{1,2}$ denoted by red dashed lines. This regime is crucial for detecting the entanglement used in RefsVitali_PRL2007 ; Palomaki_Sci2013 . in which an additional cavity has a beam-splitter-like interaction with the magnon mode for reading out the magnon states associated with the CM. The transferred entanglement can then be measured through the homodyne detection by sending a weak microwave probe. This approach requires much larger cavity leakage than the magnon dissipation, namely, $\gamma_{c}\gg\gamma_{1,2}$, so that the magnon states can remain almost unchanged when switching off the laser driving.

The time-resolved detection of the photons emitting off the cavity axis may offer an alternative scheme for entanglement measurement. The quadrature information of magnon modes can be transferred to the time-gated emitted photons, which can be homodynely detected by interfering with an extra microwave field. This quantum-light-probe scheme may take the advantage of being noninvasive detection for the entanglement measurement.

Conclusion and remarks.– In conclusion, we have proposed a protocol for entangling the magnon modes in two massive YIG spheres, through the Kerr nonlinearity that originates from the magnetocrystalline anisotropy. We shew that such nonlinearity has to be essentially included, for producing the entanglement. Our work demonstrated the stationary entanglement between two macroscopic YIG spheres driven far-from-equilibrium, within the experimentally feasible parameter regime. The amount of entanglement is quantified by the logarithmic negativity and surprisingly robust against the cavity leakage: entangled quantum state may persist with low-quality cavity giving weak magnon-cavity coupling. This may be helpful to the experimental design for the entanglement measurement.

We should note that our idea for entangling magnon modes may be potentially extensive to other complex systems, such as molecular aggregates and clusters, along with the fact of similar forms of nonlinear couplings $b^{\dagger}bq$ and $\Delta b^{\dagger}bb^{\dagger}b$. With the scaled-up parameters, the long-range entanglement in molecular aggregates would be anticipated, in that the exciton-exciton interaction is of several orders of magnitude higher than the Kerr nonlinearity resulting from the magnetocrystalline anisotropy. For instance, the two-exciton coupling in J-aggregate and light-harvesting antenna take the value of $\sim 50$cm*-1* which is $\sim 0.3\%$ of the exciton frequency. This is much stronger nonlinearity than that in YIGs with Kerr coefficient $K\sim 0.1$nHz that is $\sim 10^{-11}$ of its Kittel frequency. Recent development in both ultrafast spectroscopy and synthesis have revealed the important role of quantum coherence which may significantly modify the functions of complex molecules and may help the design of polaritonic molecular devices as well as polariton chemistry. Hence entangling the molecular aggregates may help the studies of quantum phenomena in complex molecules.

We gratefully acknowledge the support of AFOSR Award FA-9550-18-1-0141, ONR Award N00014-16-1-3054 and Robert A. Welch Foundation (Award A-1261 & A-1943-20180324). We also thank Jie Li and Tao Peng for the useful discussions.