Enhanced the Bi-(Tri-)partite entanglement between three magnons in a cavity magnomechanics system

The collective excitations of spins known as magnons in ferrimagnetic materials such as the yttrium-iron-garnet (YIG) spheres, have attracted the attention of researchers in recent years owing to their high spin density and low damping rate. The magnon mode in the YIG sphere can strongly couple to the microwave (MW) cavity photons, resulting in cavity-polaritons [1–4], which provide a potentially useful framework for magnon-based information processing, including the remote manipulation of spin current [5], the magnon dark modes and gradient memory [6] and bistability [7].

Another exciting option is to couple the magnons to the phonons [8, 9] or superconducting qubits [10–12] in order to transform magnon-cavity systems into hybrid systems. In particular, the cavity magnomechanical (CMM) system [8] consists of an MW cavity and a YIG sphere where the magnon mode in the YIG sphere couple simultaneously to the MW cavity and phonon modes through the magnetic dipole interaction and magnetostrictive force (which plays the role of the radiation pressure in optomechanics (OM) systems), respectively. Compared to the OM systems [13–19], the CMM systems benefit from more flexibility and lower loss, which offers a potential vision for implementing quantum information theory and macroscopic quantum effects. Among these effects, quantum entanglement is a primary physical resource for quantum information processing. It is an essential part of a wide variety of fields of study that include quantum dense coding [20], quantum key distribution [21], quantum teleportation [22], and quantum cryptography [23].

In the CMM system, Li et al have utilized the nonlinearity of magnetostrictive interaction to investigate various interesting quantum effects. They have demonstrated the generation of bipartite entanglements and even a genuine tripartite photon-magnon-phonon entanglement [24], which is the most significant feature of this system. Additionally, they have also demonstrated squeezing of the MW fields [25] and entangled two magnons by coupling two YIG to a single MW cavity and activating the magnomechanical interaction [26]. In the same vein, other researchers have investigated other phenomena. For instance, we may cite the generation of stationary entanglement between two MW fields by using magnetostrictive interaction [27–29] and the entanglement of two phonon modes of two YIG by utilizing the nonlinearity of the magnetostriction [30]. Interestingly enough, other nonlinearities have been used in the CMM systems [31], such as the nonlinearity of a flux-driven Josephson parametric amplifier to generate squeezed states of magnon and phonon modes [32] and entanglement of two magnon modes [33–37]. Additionally, the nonlinearity of the Kerr effect was utilized to enhance the entanglement between two magnons [38] and squeeze the state of magnon mode to enhance the entanglement [39]. Recently, the nonlinearity of the optical parametric amplifier has been utilized to improve entanglement in the CMM system [40, 41], entangle two cavity modes, and squeeze the magnon mode [42]. Moreover, a coherent feedback loop has been also employed to significantly improve entanglement in the CMM system [43].

Here, we present a model for a CMM system which involves a MW cavity field, three magnon modes and a phonon mode. Specifically, we show that by using the nonlinearity of the magnetostrictive interaction activated in one YIG sphere, it is possible to generate, both the bipartite entanglement between distinct magnon-mode subsystems and also a genuine tripartite entanglement. Furthermore, by choosing the appropriate experimental parameters in our simulations, we show that the entanglement is robust to the thermal influences. In our scheme, the entanglement results from the combination of the cavity-magnon beam-splitter interaction and the nonlinearity of magnetostrictive force, which can be enhanced by driving the YIG sphere using a laser.

This paper is organized as follows: In section 2, we introduce our system and its Hamiltonian. In section 3, we give the quantum Langevin equations (QLEs) and investigate the dynamical equation for the covariance matrix (CM), i.e. the steady-state solution of QLEs, we also give the definition of bipartite and tripartite entanglements. Section 4 is devoted to the examination of the feasibility of the system. In section 5, we present and discuss our numerical results. Finally, our conclusions are presented in section 6.

We propose in figure 1(a) hybrid CMM system which contains five modes: a MW cavity mode (a), the three magnon modes (m₁, m₂ and m₃), and a mechanical mode (b). The MW cavity field and the three magnon modes interact through the magnetic dipole interaction. Due to the high spin density of the YIG spheres, this coupling can be very strong [1, 2, 9]. In addition, the magnetostrictive force mediates the coupling between the magnon and the phonon modes of the YIG-1 sphere. In particular, the variation of magnetization inside the YIG-1 caused by the magnon excitation leads to the deformation of its geometry structure and yields the magnon1-phonon coupling [24]. Furthermore, since the dimensions of the YIG spheres are considerably smaller than the wavelength of the MW field, the influence of radiation pressure can be considered negligible [24]. The Hamiltonian of the system with respect to the magnon drive frequency ω_d reads:

Zoom In Zoom Out Reset image size

$$\begin{align} \mathcal{H}/\hbar&=\Delta _{a}a^{\dagger }a+\dfrac{\omega_{b}}{2}\left(q^{2}+p^{2}\right)\nonumber\\&\quad+\sum_{j=1,2,3}\left[\Delta_{j}m_{j}^{\dagger }m_{j}+g_{j}\left(am_{j}^{\dagger }+a^{\dagger }m_{j}\right)\right]\nonumber\\&\quad +g_{mb}m_{1}^{\dagger }{m_{1}}q +i\Omega\left(m_{1}^{\dagger} -m_{1}\right), \end{align} \tag{ 1 }$$

where $a^{\dagger }$/a and $m_{j}^{\dagger }$/m_j being the creation/annihilation operators of the cavity mode and jth magnon mode, respectively, with the corresponding frequencies ω_a and ω_j , fulfilling the commutation relation $[O,O^{\dagger}] = 1$ ($O = a,m_{j}$). p and q are dimensionless momentum and position quadratures of the mechanical mode with corresponding frequency ω_b , fulfilling the commutation relation $[q,p] = i$. In addition, $\Delta_{a} = \omega_{a}-\omega_{d}$ and $\Delta_{j} = \omega_{j}-\omega_{d}$ are the detunings, respectively, of the cavity mode and the jth magnon mode defined from the driven field. By altering the external bias magnetic field (H), the frequency of magnon modes can be adjusted over a wide range via $\omega_{j} = \gamma_{0}H$ with $\gamma_{0}/2\pi = 28$ GHz T⁻¹ being the gyromagnetic ratio. g_j denotes the coupling force between the cavity mode and the jth magnon mode, which corresponds to the strong coupling regime when g_j is greater than the decay rates of the cavity (κ_a ) and the jth magnon (κ_j ) modes, i.e. $g_{j}\gt\gt\kappa_{a}, \kappa_{j}$ [1–3]. In addition, g_mb denotes the magnomechanical coupling strength which is generally small, but it can be increased by applying a powerful MW field to excite the magnon mode, which can be realized for example by using a small loop antenna at the end of a superconducting MW line to drive the YIG sphere [7, 27, 44]. The Rabi frequency $\Omega = \frac{\sqrt{5}}{4}\gamma _{0}\sqrt{N_{t}}B_{0}$ [24] represents the coupling strength between the magnon mode and its driving magnetic field with amplitude B₀ and frequency ω_d , where $N_{t} = \rho V$ is the total spin number with the spin density $\rho = 4.22\times 10^{27}$ $\mathrm{m}^{-3}$ and V the volume of the sphere [24]. It is worth noting here that the bosonic operators $m_{j}^{\dagger}$ and m_j are derived under the assumption of low-lying excitations, i.e. $2N_{t}s\gt\gt\langle m_{j}^{\dagger}m_{j} \rangle$, where $s = \dfrac{5}{2}$ is the spin number of the ground-state Fe³⁺ ion in the YIG sphere [24].

In this section, we derive the dissipative dynamics of our system described by the following set of QLEs:

$$\begin{align} \dot{a} =-(i\Delta _{a}+\kappa _{a})a-ig_{1}m_{1}-ig_{2}m_{2}-ig_{3}m_{3}+ \sqrt{2\kappa _{a}}a^{\mathrm{in}}, \end{align} \tag{ 2 }$$

$$\begin{align} \kern-9.5pt\dot{m}_{1} =-(i\Delta _{1}+\kappa _{1})m_{1}-ig_{1}a-ig_{mb}m_{1}q+\Omega + \sqrt{2\kappa _{1}}m_{1}^{\mathrm{in}}, \nonumber\\ \end{align} \tag{ 3 }$$

$$\begin{align} \kern-72.1pt\dot{m}_{2} =-(i\Delta _{2}+\kappa _{2})m_{2}-ig_{2}a+\sqrt{2\kappa _{2}} m_{2}^{\mathrm{in}}, \end{align} \tag{ 4 }$$

$$\begin{align} \dot{m}_{3} =-(i\Delta _{3}+\kappa _{3})m_{3}-ig_{3}a+\sqrt{2\kappa _{3}} m_{3}^{\mathrm{in}}, \phantom{0000000000} \end{align} \tag{ 5 }$$

$$\begin{align} \dot{q} =\omega _{b}p, \phantom{0000000000000000000000000000000000} \end{align} \tag{ 6 }$$

$$\begin{align} \kern-60.2pt\dot{p} =-\omega _{b}q-\kappa _{b}p-g_{mb}m_{1}^{\dagger }m_{1}+\xi ,\phantom{000000} \end{align} \tag{ 7 }$$

where κ_b is the mechanical dissipation rate, while $a^{\mathrm{in}}$ and $m_j^{\mathrm{in}}$ denote the input noise operators for the cavity and the jth magnon modes, respectively. These operators have a mean of zero and are associated with the following correlation functions [45]:

$$\begin{align} \;\;\left\langle a^{in\dagger }(t)a^{\mathrm{in}}(t^{^{\prime} });a^{\mathrm{in}}(t)a^{in\dagger }(t^{^{\prime} }) \right\rangle =[ n_{a}(\omega _{a});n_{a}(\omega _{a})+1]\delta (t-t^{^{\prime} }), \end{align} \tag{ 8 }$$

$$\begin{align} \left\langle m_{j}^{in\dagger }(t)m_{j}^{\mathrm{in}}(t^{^{\prime} });m_{j}^{\mathrm{in}}(t)m_{j}^{in\dagger}(t^{^{\prime} }) \right\rangle =[ n_{j}(\omega _{j});n_{j}(\omega _{j})+1]\delta (t-t^{^{\prime} }),\quad \end{align} \tag{ 9 }$$

$$\begin{align} \phantom{0}\left\langle \xi (t)\xi (t^{^{\prime} })+\xi (t^{^{\prime} })\xi (t)\right\rangle /2 &\simeq \kappa _{b}[2n_{b}(\omega _{b})+1]\delta (t-t^{^{\prime} }),\; \end{align} \tag{ 10 }$$

where the correlation function of the Langevin force operator ξ has been taken into account for the mechanical Brownian motion under the Markovian approximation, which is valid for a large mechanical quality factor, i.e. $Q = \dfrac{\omega_{b}}{\kappa_{b}}\gt\gt1$ [46, 47]. $n_{O}(\omega _{O}) = [\exp (\frac{\hbar \omega _{O}}{k_\mathrm{B}T})-1]^{-1}$ ($O = a,j,b$) being the equilibrium mean thermal photon, jth magnon, and phonon numbers, with $k_\mathrm{B}$ is the Boltzmann constant and T is the environmental temperature.

Since the first magnon mode (m₁) is strongly driven, its steady-state amplitude is much larger than one, i.e. $\left|\left\langle m_{1}\right\rangle \right| \gt\gt1$. This large amplitude of the first magnon mode leads to significant amplitudes of the cavity, magnon 2, and magnon 3 modes, respectively, due to the beam-splitter interaction between the magnons and the cavity mode, which implies that the amplitudes of $\left|\left\langle a\right\rangle \right|, \left|\left\langle m_{k}\right\rangle \right| \gt\gt1$. As a result, the QLEs (2)–(7) can be linearized using the standard linearization technique. This involves expressing each Heisenberg operator as a sum of its classical average and a fluctuation operator with zero mean, i.e. $O = \left\langle O\right\rangle+ \delta O$ ($O = a,m_{j},q,p$), and neglecting higher-order fluctuations. After linearization, two types of linearized QLEs are obtained for the classical averages and fluctuation operators. The solutions for the first type of equations are:

$$\begin{align} \langle a\rangle &\simeq -\frac{i\Omega g_{1}\Delta _{2}\Delta _{3}}{g_{3}^{2}\tilde{\Delta} _{1}\Delta _{2}+g_{2}^{2}\tilde{\Delta}_{1}\Delta _{3}+g_{1}^{2}\Delta _{2}\Delta _{3}-g_{1}^{2}\Delta _{a}\tilde{\Delta}_{1}\Delta _{2}\Delta _{3}}, \end{align} \tag{ 11 }$$

$$\begin{align} \langle m_{1}\rangle &\simeq -\frac{i\Omega \left[ g_{3}^{2}\Delta _{2}+\Delta _{3}(g_{2}^{2}-\Delta _{a}\Delta _{2})\right] }{ g_{3}^{2}\tilde{\Delta}_{1}\Delta _{2}+g_{2}^{2}\tilde{\Delta}_{1}\Delta _{3}+g_{1}^{2}\Delta _{2}\Delta _{3}-g_{1}^{2}\Delta _{a}\tilde{\Delta}_{1}\Delta _{2}\Delta _{3}}, \end{align} \tag{ 12 }$$

$$\begin{align} \langle m_{2}\rangle &\simeq -\frac{i\Omega g_{1}g_{2}\Delta _{2}}{g_{3}^{2}\tilde{\Delta}_{1}\Delta _{2}+g_{2}^{2}\tilde{\Delta}_{1}\Delta _{3}+g_{1}^{2}\Delta _{2}\Delta _{3}-g_{1}^{2}\Delta _{a}\tilde{\Delta}_{1}\Delta _{2}\Delta _{3}}, \end{align} \tag{ 13 }$$

$$\begin{align} \langle m_{3}\rangle &\simeq -\frac{i\Omega g_{1}g_{3}\Delta _{2}}{g_{3}^{2}\tilde{\Delta}_{1}\Delta _{2}+g_{2}^{2}\tilde{\Delta}_{1}\Delta _{3}+g_{1}^{2}\Delta _{2}\Delta _{3}-g_{1}^{2}\Delta _{a}\tilde{\Delta}_{1}\Delta _{2}\Delta _{3}}, \end{align} \tag{ 14 }$$

$$\begin{align} \langle p\rangle =0, \langle q\rangle = - \frac{g_{mb}}{\omega _{b}}\left\vert \langle m_{1}\rangle \right\vert^{2},\phantom{000000000000000000} \end{align} \tag{ 15 }$$

where the above expressions are reached within the condition that $\vert \Delta_{a} \vert, \vert \tilde{\Delta}_{1} \vert,\vert \Delta_{k} \vert \gt\gt \kappa_{a},\kappa_{1},\kappa_{k}$ ($k = 2,3$) and $\tilde{\Delta}_{1} = \Delta_{1}+g_{mb} \langle q\rangle$ represents the effective detuning of the first magnon mode, which takes into account the frequency shift brought on by the interaction between magnon1 and mechanical mode.

Alternatively, the linearized QLEs can be represented in quadrature form as:

$$\begin{align} \delta \dot{X} &=-\kappa _{a}\delta X+\Delta _{a}\delta Y+g_{1}\delta y_{1}+g_{2}\delta y_{2}+g_{3}\delta y_{3}+\sqrt{2\kappa _{_{a}}}X^{\mathrm{in}},\; \end{align} \tag{ 16 }$$

$$\begin{align} \delta \dot{Y} =-\kappa _{a}\delta Y-\Delta _{a}\delta X-g_{1}\delta x_{1}-g_{2}\delta x_{2}-g_{3}\delta x_{3}+\sqrt{2\kappa _{a}}Y^{\,\mathrm{in}}, \end{align} \tag{ 17 }$$

$$\begin{align} \delta \dot{x}_{1} =-\kappa _{1}\delta x_{1}+\tilde{\Delta}_{1}\delta y_{1}+g_{1}\delta Y-G_{mb}\delta q+\sqrt{2\kappa _{1}}x_{1}^{\mathrm{in}},\quad \end{align} \tag{ 18 }$$

$$\begin{align} \delta \dot{y}_{1} =-\kappa _{1}\delta y_{1}-\tilde{\Delta}_{1}\delta x_{1}-g_{1}\delta X+\sqrt{2\kappa _{1}}y_{1}^{\mathrm{in}}, \qquad\qquad\; \end{align} \tag{ 19 }$$

$$\begin{align} \delta \dot{x}_{k} =-\kappa _{k}\delta x_{k}+\Delta _{k}\delta y_{k}+g_{k}\delta Y+\sqrt{2\kappa _{k}}x_{k}^{\mathrm{in}},\text{}k=2,3\quad\quad \end{align} \tag{ 20 }$$

$$\begin{align} \delta \dot{y}_{k} =-\kappa _{k}\delta y_{k}-\Delta _{k}\delta x_{k}-g_{k}\delta X+\sqrt{2\kappa _{k}}y_{k}^{\mathrm{in}},\text{}k=2,3\quad\quad \end{align} \tag{ 21 }$$

$$\begin{align} \delta \dot{q} =\omega _{b}\delta p,\phantom{0000000000000000000000000000000} \end{align} \tag{ 22 }$$

$$\begin{align} \delta \dot{p} =-\kappa _{b}\delta p-\omega _{b}\delta q+G_{mb}\delta y_{1}+\xi ,\qquad\qquad\qquad\quad \end{align} \tag{ 23 }$$

in which we have characterized the quadrature fluctuations defined as: $\delta X = \frac{1}{\sqrt{2}}\left( \delta a^{\dagger }+\delta a\right), \delta Y = \frac{i}{\sqrt{2}}\left( \delta a^{\dagger }-\delta a\right),\delta x_{j} = \frac{1}{\sqrt{2}} \left( \delta m_{j}^{\dagger }+\delta m_{j}\right), \delta y_{j} = \frac{i}{\sqrt{2}}\left( \delta m_{j}^{\dagger }-\delta m_{j}\right)$. $G_{mb} = i\sqrt{2}g_{mb}\left\langle m \right\rangle$ is the effective magnomechanical coupling strength. In addition, the aforementioned linearized QLEs (16)–(23) can be expressed in the following compact matrix form:

$$\begin{align} \dot{\mathcal{U}}(t)=\mathcal{A}\mathcal{U}(t)+\mathscr{N}(t), \end{align} \tag{ 24 }$$

with $\mathcal{U}^\top(t)\! =\! [ \delta X(t), \delta Y(t),\delta x_{1}(t),\delta y_{1}(t) ,\delta x_{2}(t), \delta y_{2}(t), \delta x_{3}(t),$ $\delta y_{3}(t), \delta q(t), \delta p(t)]$ and $\mathscr{N}^\top(t) = [ \sqrt{2\kappa _{a}} X^{\mathrm{in}}(t), \sqrt{2\kappa _{a}} Y^{\,\mathrm{in}}(t),$ $\sqrt{2\kappa _{1}}x_{1}^{\mathrm{in}}(t),$ $\sqrt{2\kappa _{1}}y_{1}^{\mathrm{in}}(t),$ $\sqrt{2\kappa _{2}}x_{2}^{\mathrm{in}}(t),$ $\sqrt{2\kappa _{2}}y_{2}^{\mathrm{in}}(t),$ $\sqrt{2\kappa _{3}}x_{3}^{\mathrm{in}}(t)$,

$\sqrt{2\kappa _{3}}y_{3}^{\mathrm{in}}(t), 0, \xi(t)]$ are the vectors for quantum fluctuations and noise, respectively, and the matrix $\mathcal{A}$ reads:

$$\begin{align} \mathcal{A}= \begin{pmatrix} -\kappa _{a} & \Delta _{a} & 0 & g_{1} & 0 & g_{2} & 0 & g_{3} & 0 & 0 \\ -\Delta _{a} & -\kappa _{a} & -g_{1} & 0 & -g_{2} & 0 & -g_{3} & 0 & 0 & 0 \\ 0 & g_{1} & -\kappa _{1} & \tilde{\Delta} _{1} & 0 & 0 & 0 & 0 & -G_{mb} & 0 \\ -g_{1} & 0 & -\tilde{\Delta} _{1} & -\kappa _{1} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & g_{2} & 0 & 0 & -\kappa _{2} & \Delta _{2} & 0 & 0 & 0 & 0 \\ -g_{2} & 0 & 0 & 0 & -\Delta _{2} & -\kappa _{2} & 0 & 0 & 0 & 0 \\ 0 & g_{3} & 0 & 0 & 0 & 0 & -\kappa _{3} & \Delta _{3} & 0 & 0 \\ -g_{3} & 0 & 0 & 0 & 0 & 0 & -\Delta _{3} & -\kappa _{3} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \omega _{b} \\ 0 & 0 & 0 & G_{mb} & 0 & 0 & 0 & 0 & -\omega _{b} & -\kappa _{b} \end{pmatrix}. \end{align} \tag{ 25 }$$

Since we are working with linearized versions of the QLEs and all noises are Gaussian, the quantum fluctuation steady-state of the system is a continuous variable Gaussian state with five modes. This state can be described by a $10\times10$ CM whose elements are defined by $\mathcal{V}_{kl} = \left\langle \mathcal{U}_{k}(\infty)\mathcal{U}_{l}(\infty)+\mathcal{U}_{l}(\infty)\mathcal{U}_{k}(\infty)\right\rangle /2$ for $k,l = 1,2,{\ldots},10$.

To derive the steady-state CM $\mathcal{V}$, we can solve the following Lyapunov equation, as guaranteed by the negativity of the real parts of the eigenvalues of the matrix $\mathcal{A}$ [48]:

$$\begin{align} \mathcal{A}\mathcal{V}+\mathcal{V}\mathcal{A}^\top=-\mathcal{R}, \end{align} \tag{ 26 }$$

where $\mathcal{R}$ is the diffusion matrix of elements: $\mathcal{R}_{kl}\delta(t-t^{^{\prime}}) = (\left\langle \mathscr{N}_{k}(t)\mathscr{N}_{l}(t^{^{\prime} })+\mathscr{N}_{l}(t^{^{\prime}})\mathscr{N}_{k}(t)\right\rangle )/2$ ($k,l = 1,2,{\ldots},10$) and given by: $\mathcal{R} = {\mathrm{diag}}(\kappa_{\mathrm{a}}(2\mathrm{n}_{\mathrm{a}}+1),\kappa_{\mathrm{a}}(2\mathrm{n}_{\mathrm{a}}+1), \kappa_{1}(2\mathrm{n}_{1}+1),\kappa_{1}(2\mathrm{n}_{1}+1),$ $\kappa_{2}(2\mathrm{n}_{2}+1),$ $\kappa_{2}(2\mathrm{n}_{2}+1),$ $\kappa_{3}(2\mathrm{n}_{3}+1),$ $\kappa_{3}(2\mathrm{n}_{3}+1),0,\kappa_{\mathrm{b}}(2\mathrm{n}_{\mathrm{b}}+1))$.

After obtaining the CM $\mathcal{V}$, the bipartite entanglement between the different magnons and cavity mode, can be quantified using the logarithmic negativity E_N [49], which is defined as:

$$\begin{align} E_{N}=\max\left[ 0,-\mathrm{ln} (2\tilde{\eta}^{-})\right], \end{align} \tag{ 27 }$$

where $\tilde{\eta}^{-} = \min\mathrm{eig}|i\mathcal{J}_{2}\tilde{\mathcal{V}}_{4}|$ is the smallest symplectic eigenvalue of the partial transpose CM $\tilde{\mathcal{V}}_{4} = \mathcal{P}_{1/2}\mathcal{V}_{4} \mathcal{P}_{1/2}$. $\mathcal{J}_{2} = \bigoplus_{j = 1}^{2}i\sigma_{2}$ (with σ₂ is the y-Pauli matrix) is the symplectic matrix, $\mathcal{P}_{1/2} = \mathrm{diag}(1,-1,1,1)$ is the matrix that performs the partial transposition at the level of the CMs and $\mathcal{V}_{4}$ is the $4 \times4$ CM of the two interested subsystems, which derived by deleting the uninteresting mode's rows and columns in $\mathcal{V}$. In this CMM system, the presence of bipartite entanglement is reflected by the nonzero of the logarithmic negativity, i.e. $E_{N}\gt0$.

Next, we investigate the tripartite entanglement $m_{1} - m_{2} - m_{3}$ ($a - m_{1} - m_{2}$) as measured by the residual contangle, which is given by [24, 38, 49, 50]:

$$\begin{align} \mathcal{R}_{i/jk }=\mathcal{C}_{i/jk}-\mathcal{C}_{i/j}-\mathcal{C}_{i/k}, \end{align} \tag{ 28 }$$

where the indices i, j, and k refer to three modes a, m₁, and m₂, or to m₁, m₂, and m₃, respectively. $\mathcal{C}_{u/v}$ represents the contangle between subsystems u and v. In this context, subsystem v may consist of one or two modes. The contangle is a proper entanglement monotone defined as the squared logarithmic negativity [49, 50], i.e. $\mathcal{C}_{u/v} = E^{2}_{u/v}$. The logarithmic negativity of one mode with respect to that of two modes $E_{i/jk}$ is defined as follows:

$$\begin{align} E_{i/jk}=\max\left[ 0,-\mathrm{ln } 2\tilde{\eta}^{-}_{i/jk}\right], \end{align} \tag{ 29 }$$

similar to the logarithmic negativity for bipartite entanglement (27), $\tilde{\eta}_{i/jk}^{-} = \min \mathrm{eig}|i\mathcal{J}_{3}\tilde{\mathcal{V}}_{6}|$ is the smallest symplectic eigenvalue of the $6\times 6$ CM $\tilde{\mathcal{V}}_{6} = \mathcal{P}_{i/jk}\mathcal{V}_{6} \mathcal{P}_{i/jk}$ and $\mathcal{J}_{3} = \bigoplus_{j = 1}^{3}i\sigma_{2}$ being the symplectic matrix. Here $\mathcal{V}_{6}$ is the $6\times 6$ CM describing the three interesting modes which can be obtained by deleting the uninteresting mode's rows and columns in $\mathcal{V}$. $\mathcal{P}_{i/jk} = \sigma _{3}\oplus \unicode{x1D7D9}\oplus \unicode{x1D7D9}$, $\mathcal{P}_{j/ik} = \unicode{x1D7D9}\oplus\sigma _{3} \oplus \unicode{x1D7D9}$ and $\mathcal{P}_{k/ij} = \unicode{x1D7D9}\oplus\unicode{x1D7D9} \oplus \sigma _{3}$ are partial transposition matrices at the level of the CM $\mathcal{V}_{6}$. Furthermore, the residue contangle fulfills the monogamy of quantum entanglement, that is, $\mathcal{R}_{i/jk }\geqslant0$, i.e. $\mathcal{C}_{i/jk}\geqslant\mathcal{C}_{i/j}+\mathcal{C}_{i/k}$, which is similar to the monogamy inequality for a three-qubit cases established by Coffman et al in [51]. A bona fide minimum residual cotangle $\mathcal{R}_\mathrm{min}$ provides a genuine measure of tripartite entanglement for Gaussian states, which is given by:

$$\begin{align} \mathcal{R}_{\min}=\min[\mathcal{R}_{i/jk},\mathcal{R}_{j/ik},\mathcal{R}_{k/ij}]. \end{align} \tag{ 30 }$$

This relation ensures the invariance of tripartite entanglement under any permutation of the three modes. A nonzero minimal residual contangle $\mathcal{R}_{\min}\gt0$ indicates that the three modes are genuinely entangled.

Let's evaluate the feasibility of our proposed scheme by utilizing experimentally feasible parameters taken from [24]. The parameters are as follows: $\kappa_{a}/2\pi = \kappa_{1}/2\pi = \kappa_{2}/2\pi = \kappa_{3}/2\pi = 1$ MHz, $\kappa_{b}/2\pi = 100$ Hz, $\omega_{a}/2\pi = 10$ GHz, $\omega_{b}/2\pi = 10$ MHz, T = 10 mK, $G_{mb}/2\pi = 4.8$ MHz and $g_{1}/2\pi = 4.8$ MHz.

On one hand, we made the assumption of a low-lying excitation for which the numerical results are applicable, i.e. when the magnon excitation numbers satisfy $\langle m_{j}^{\dagger}m_{j}\rangle \lt\lt2N_{t}s = 5N_t$ [24]. According to references [24, 26, 30], in our simulation, we used YIG spheres with a diameter of 250 $\mu m$, containing $N_{t} \simeq 3.5\times 10^{16}$ spins, and $G_{mb}/2\pi = 4.8\times 10^{6}$ Hz, corresponding to a Rabi frequency $\Omega\simeq 7.12\times 10^{14}$ Hz. Additionally, $\left|\langle m_{1}\rangle \right| \simeq 1.1\times 10^{7}$, so $\langle m_{1}^{\dagger}m_{1}\rangle \simeq 1.3\times 10^{14}\lt\lt2N_ts = 5N_t = 1.7\times 10^{17}$, which is holds. Similarly, $\langle m_{2}^{\dagger}m_{2}\rangle \lt\lt2N_ts$ and $\langle m_{3}^{\dagger}m_{3}\rangle \lt\lt2Ns$ are also satisfied for the parameters used in this manuscript because $|\langle m_{2}^{\dagger}m_{2}\rangle |, | \langle m_{3}^{\dagger}m_{3}\rangle | \lt\lt|\langle m_{1}^{\dagger}m_{1}\rangle |$.

On the other hand, if a strong pump is utilized to drive the magnon mode, it can introduce the Kerr nonlinear term $K m^{\dagger}m m^{\dagger}m$ into the Hamiltonian [7, 52], where the value of the Kerr coefficient K is determined by the volume of the YIG sphere and is inversely proportional to it. For a YIG sphere with a diameter of 250 $\mu\mathrm{m}$ [24], the Kerr coefficient K is estimated to be approximately $2\pi\times6.4 \times 10^{-9}$ Hz [24]. The Kerr effect can be considered negligible only if the condition $K|\langle m_{1}\rangle|^{3}\lt\lt\Omega$ is well satisfied. The parameters used in our simulations lead to $K | \langle m_1\rangle|^{3}\simeq 5.8\times 10^{13}$ Hz $\lt\lt$ $\Omega\simeq 7.12\times 10^{14}$ Hz. Therefore, we can conclude that the Kerr nonlinearity can be neglected in our scheme.

In this section, we will explore the results of the study regarding the entanglement between different magnons and cavity mode, including bipartite and tripartite entanglement.

In figures 2(a) and (b), we present the bipartite entanglements, namely (a) cavity-magnon1 $E_{am_1}$ and (b) magnon1-phonon $E_{m_1b}$, as a functions of the normalized detunings $\Delta_{a}/\omega_{b}$ and $\tilde{\Delta}_{1}/\omega_{b}$. It is important to mention that all results in this paper are in steady-state, which has been verified by ensuring that the real parts of the eigenvalues of the matrix $\mathcal{A}$ are negative (see appendix). It can be seen from figure 2(a) that there is a specific set of detunings, centered around $\tilde{\Delta}_{1}\simeq \omega_{b}$ and $\Delta_{a}\simeq- \omega_{b}$, where the magnon1 and cavity modes are in resonance with the two mechanical sidebands. This detunings implies strong entanglement between the magnon1 and cavity modes. Furthermore, applying a strong red-detuned driving field to the magnon1 mode results in the cooling of the hot mechanical mode b (which is a necessary step in preparing quantum states), as shown in [24]. At the same time, this activates the magnomechanical parametric down-conversion interaction and generates significant magnomechanical entanglement $E_{m_1b}$ due to the strong coupling G_mb (see figure 2(b)) [24]. This results in partial transfer of entanglement to the cavity-magnon1 subsystem, meaning that both the cavity and magnon1 mode become entangled, as evidenced by a positive value of $E_{am_{1}} \gt 0$, and the same effect is observed for the cavity-phonon subsystem [24]. This entanglement transfer is more pronounced when the cavity detuning is approximately $\Delta_{a}\simeq-\omega_{b}$ (see figure 2(a)). On the other hand, we plot in figure 2(c) $E_{am_1}$ versus $\Delta_{a}/\omega_{b}$ and $G_{mb}/g_{1}$. As can be seen from figure 2(c), when $G_{mb} = 0$, we note that the bipartite entanglement $E_{am_{1}} = 0$, which is natural due to the fact that the linear cavity-magnon beam splitter interaction generally generates zero entanglement, and only the nonlinearity of the magnetostrictive force is responsible for generating the squeezing and entanglement in the CMM system [24, 32]. We also note that strong entanglement is generated when the two couplings are of the same magnitude, i.e. $G_{mb}\simeq g_{1}$.

Zoom In Zoom Out Reset image size

Figure 3 presents the bipartite cavity-magnon2 entanglement ($E_{am_2}$) as a function of (a) $\Delta_{a}/\omega_{b}$ and $\tilde{\Delta}_{1}/\omega_{b}$, (b) $\Delta_{a}/\omega_{b}$ and $\Delta_{2}/\omega_{b}$, and (c) $g_{2}/g_{1}$ and $g_{3}/g_{1}$. Figures 3(a) and (b) show that there are optimal detunings for the cavity and magnon modes that lead to strong entanglement between the magnon2 and cavity modes. These detunings are given as $\Delta_{a}\simeq0$ (i.e. when the cavity mode resonates with the driving field), $\tilde{\Delta}_{1}\simeq0.8\omega_{b}$, and $\Delta_{2}\simeq-0.7\omega_{b}$. This entanglement is achieved through a combination of nonlinear magnetostrictive interaction and linear beam splitter interaction, as well as carefully chosen detunings for each mode. The nonlinearity of the magnetostrictive interaction is essential for creating entangled states that cannot be achieved using linear interaction alone (see figure 2), while the linear interaction enhances the entanglement and allows for efficient information transfer between the cavity and jth magnon subsystems. Furthermore, figure 3(c) shows that there are also optimal couplings that lead to significant entanglement between the magnon2 and cavity modes when $g_{2} = 0.8g_{1}$ and $g_{3} = 0.3g_{1}$. It is worth mentioning that similar optimal conditions for strong entanglement between the cavity and magnon3 modes can be also obtained by interchanging the roles of m₂ and m₃ in the analysis. However, to keep the discussion concise, we have not explored this scenario in detail here.

Zoom In Zoom Out Reset image size

In order to examine the entanglement between the three magnon modes, we have presented bipartite entanglement $E_{m_{1}m_{2}}$ in figure 4 and analyzed it against (a) $\Delta_{a}/\omega_{b}$ and $\tilde{\Delta}_{1}/\omega_{b}$, (b) $\Delta_{a}/\omega_{b}$ and $\Delta_{2}/\omega_{b}$ and (c) $\Delta_{a}/\omega_{b}$ and $g_{2}/g_{1}$. Based on figures 4(a) and (b), we can conclude that the maximum entanglement $E_{m_{1}m_{2}}$ can be achieved when the cavity and magnon2 modes are approximately resonant with the red sideband, i.e. $\Delta_{a} = \Delta_{m_2}\simeq-\omega_{b}$, while the magnon m₃ and m₁ modes are approximately resonant with the blue sideband of the mechanical mode, i.e. $\Delta_{m_1} = \Delta_{m_3}\simeq\omega_{b}$. This entanglement is possible due to the magnon-cavity interaction between the three magnon and cavity modes, which is partially transferred from the cavity-magnon1 subsystem to the cavity-magnon2 and cavity-magnon3 subsystems. It is worth mentioning that a similar scenario can be seen in achieving strong entanglement between two MW cavities using the magnetostrictive interaction as reported in [27], where the MW 1 must be resonant with the red sideband, and the MW 2 must be resonant with the blue sideband, and vice versa. Additionally, figure 4(c) demonstrates that there is an optimal coupling g₂ for achieving strong entanglement when g₂ is of the same order as g₃, i.e. $g_{2} = g_{3}\simeq 0.5g_{1}$. It is important to note that the same situation can occur for the bipartite entanglement $E_{m_1m_3}$ when interchanging the roles of m₂ and m₃, but we have not discussed it here for the sake of brevity.

Zoom In Zoom Out Reset image size

In figure 5, we display the bipartite entanglement $E_{m_{2}m_{3}}$ between two indirectly coupled magnon modes as a function of three different parameters: (a) the ratios of $\Delta_{a}/\omega_{b}$ and $\tilde{\Delta}_{1}/\omega_{b}$, (b) the ratios of $\Delta_{2}/\omega_{b}$ and $\Delta_{3}/\omega_{b}$, and (c) the ratio of coupling rates $g_{2}/g_{1}$ and $g_{3}/g_{1}$. Figures 5(a) and (b) show that strong entanglement occurs under the same detuning conditions as those illustrated in figure 4, precisely when $\Delta_{a} = \Delta_{2}\simeq -\omega_{b}$ and $\tilde{\Delta}_{1} = \Delta_{3}\simeq \omega_{b}$ for $g_{2} = 0.3g_{1}$ and $g_{3} = g_{1}$. We note here that when $g_{2} = g_{1}$ and $g_{3} = 0.3g_{1}$, the optimal detunings are $\Delta_{a} = \Delta_{3}\simeq -\omega_{b}$ and $\tilde{\Delta}_{1} = \Delta_{2}\simeq \omega_{b}$. Finally, figure 5(c) investigates the optimal values of the coupling rates g₂ and g₃, revealing that the magnon mode resonating with the red (blue) mechanical sideband should have a reduced (increased) coupling rate compared to the other magnon mode. For instance, $g_2 = g_1$ and $g_3 = 0.3g_1$ are optimal to achieve strong entanglement between the two magnons as shown in figure 5(c).

Zoom In Zoom Out Reset image size

To test the robustness of entanglement between bipartite subsystems against environmental temperature, we plotted three measures of entanglement, namely $E_{am_2}$, $E_{m_1m_2}$, and $E_{m_2m_3}$, against T in figure 6. We used the optimal parameters for each subsystem from figures 3–5 and considered two cases: $\kappa_{1} = \kappa_{2} = \kappa_{a}$ (dashed line) and $\kappa_{1} = \kappa_{2} = 0.5\kappa_{a}$ (solid line). Figure 6 reveals that all the subsystems are most robust when the damping rates of the magnon1 and magnon2 modes are low. Additionally, $E_{m_1m_2}$ is more temperature-resistant than $E_{am_2}$ and $E_{m_2m_3}$. Specifically, for $\kappa_{1} = \kappa_{2} = 0.5\kappa_{a}$ ($\kappa_{1} = \kappa_{2} = \kappa_{a}$), $E_{am_2} = 0$ at $T\simeq167$ mK ($T\simeq159$ mK), $E_{m_1m_2} = 0$ at $T\simeq237$ mK ($T\simeq209$ mK), and $E_{m_2m_3} = 0$ at $T\simeq169$ mK ($T\simeq150$ mK). Thus, the damping rate of the magnon1 and magnon2 modes can also be used to control entanglement in this system. Furthermore, the entanglement decreases with increasing temperature due to the decoherence effect. Note that the highest temperature observed here is higher than the temperature obtained in [26] and lower than the temperature of 250 mK recently observed in [28] using two indirectly coupled MW fields.

Zoom In Zoom Out Reset image size

In addition to the generation of bipartite entanglement between the different magnon modes, it is also possible to generate genuine tripartite entanglement in the system by exploiting the nonlinearity of the magnetostrictive interaction, as illustrated by the nonzero minimum residual contangle $\mathcal{R}_\mathrm{min}$ in figures 7 and 8. In figure 7 we present the tripartite cavity-magnon1-magnon2 entanglement as a function of $\Delta_a/\omega_{b}$ and temperature T for the optimal parameters obtained in figure 3. As can be seen from figure 7, the system exhibits a genuine tripartite entanglement against the detuning $\Delta_{a}$ and the tripartite entanglement survives up to the temperature of 160 mK. It is worth mentioning again that similar results for tripartite cavity-magnon1-magnon3 entanglement can be obtained by interchanging the roles of m₂ and m₃ in the analysis. However, to keep the discussion concise, we have not explored this scenario in detail here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the other hand we plot in figure 8 the tripartite magnon1-magnon2-magnon3 entanglement as a function of $\Delta_a/\omega_{b}$ [panels (a) and (b)] and temperature T [panels (c) and (d)]. We used the optimal parameters obtained from figures 4, and 5 in [panels (a) and (c)] and [panels (b) and (d)], respectively. We can see from figures 8(a) and (b) that in both cases, the system exhibits genuine tripartite entanglement against the detuning $\Delta_{a}$. Furthermore, the tripartite entanglement shows robustness against the temperature T: the best robustness highlighted when we use the optimum parameters used to achieve strong entanglement between magnon2 and magnon3. In fact, $\mathcal{R}_\mathrm{min} = 0$ for $T\simeq 112$ mK [panel (c)] and $T\simeq 156$ mK [panel (d)].

Before concluding this section, and following references [24, 26, 30], we will discuss the experimental implementation. The verification of the generated bipartite or tripartite entanglement can be achieved through measurements of the corresponding CMs [54, 55]. The cavity field quadratures can be directly measured by homodyning the cavity output, and the state of the three magnons can be indirectly measured by homodyning the cavity output of a probe field introduced into the system [24, 26, 30]. In addition, an extra optical cavity can be employed to couple the YIG-1 sphere, enabling the measurement of mechanical quadratures [24, 30, 56].

We studied a CMM system consisting of an MW cavity field, three magnons and a phonon mode, where the three magnon modes in the YIG spheres interact with the MW cavity field through the magnetic dipole interaction, and one of them additionally interacts with the phonon mode of the YIG sphere through the magnetostrictive interaction. By selecting experimentally reachable parameters and fully exploiting the aforementioned nonlinearity of the magnetostrictive force and the linear cavity-magnon (beam-splitter) interaction, we have shown that the three magnons become mutually entangled, and the steady-state of the system displays a genuine tripartite entanglement, where the entanglement is resistant to the environmental temperature. Furthermore, our findings show that to have a strong bipartite magnon-magnon entanglement, the magnon that resonates with the red (blue) mechanical sideband must have a reduced (increased) coupling rate compared to the other magnon mode. Due to the macroscopic scale of the YIG spheres ($250\,\mu\mathrm{m}$-diameter), our system exhibits macroscopic quantum states that are useful for investigating quantum/classical transitions, as well as testing decoherence theory. Therefore, research involving the creation of macroscopic entangled states may benefit from this work.

All data that support the findings of this study are included within the article (and any supplementary files).

The stability of strongly magnomechanical systems is a challenging issue that requires careful analysis. Here, we will focus on the stability of our particular system. Specifically, for our system to be stable, all the eigenvalues of the drift matrix $\mathcal{A}$ should have negative real parts. If the real part of any eigenvalue changes sign, the system becomes unstable. We can employ the Routh-Hurwitz criterion [53] to gain a more comprehensive comprehension of the range of parameters where stability takes place. As the process of calculating eigenvalues through analytical methods is complex, we have chosen to utilize a numerical method to address this issue.

In figure 9, we present the maximum real parts of the eigenvalues of $\mathcal{A}$, denoted as $\lambda_{\max}$, plotted against the normalized detunings. We define $\lambda_{\max}$ as $\lambda_{\max} = \max[\text{real}(\lambda_{\mathcal{A}})]$ with $|\mathcal{A}-\lambda_{\mathcal{A}}\unicode{x1D7D9}| = 0$. It is important to note that if $\lambda_{\max}\geqslant0$, the system will become unstable. However, according to figure 9, it is clear that the system remains stable for the given parameters as the maximum eigenvalue $\lambda_{\max}$ remains negative.

Zoom In Zoom Out Reset image size

Thus, the stability conditions are satisfied by all numerical parameters used in this manuscript, leading to the conclusion that the selected working regime is stable.