††thanks<EMAIL_ADDRESS>
# Magnon squeezing via reservoir-engineered optomagnomechanics

Zhi-Yuan Fan  Huai-Bing Zhu  Hao-Tian Li  Jie Li  Zhejiang Key Laboratory of
Micro-nano Quantum Chips and Quantum Control, School of Physics, and State Key
Laboratory for Extreme Photonics and Instrumentation, Zhejiang University,
Hangzhou 310027, China

###### Abstract

We show how to prepare magnonic squeezed states in an optomagnomechanical
system, in which magnetostriction induced mechanical displacement couples to
an optical cavity via radiation pressure. We discuss two scenarios depending
on whether the magnomechanical coupling is linear or dispersive. We show that
in both cases the strong mechanical squeezing obtained via two-tone driving of
the optical cavity can be efficiently transferred to the magnon mode. In the
linear coupling case, stationary magnon squeezing is achieved; while in the
dispersive coupling case, a transient magnonic squeezed state is prepared in a
two-step protocol. The proposed magnonic squeezed states find promising
applications in quantum information processing and quantum sensing using
magnons.

##  I Introduction

The past decade has witnessed the formation and fast development of cavity
magnonics, review1 ; review2 ; review3 ; review4  which studies the
interaction between microwave cavity photons with collective spin excitations
(magnons) in magnetic materials, such as yttrium iron garnet (YIG). One of the
most prominent advantages of the magnonic system is that it exhibits an
excellent ability to coherently couple with a variety of physical systems,
including microwave photons, Huebl2013 ; Tabuchi2014 ; Zhang2014  optical
photons, Osada2016 ; Tang2016 ; Haigh2016  vibration phonons, Zhang2016 ;
Li2018 ; Potts2021 ; Shen2022  and superconducting qubits, Tabuchi2015 ;
DLQ2017 ; DLQ2020 ; Xu2023  etc. The magnonic system therefore has great
potential for building a hybrid quantum system with multiple components of
different nature, which finds broad applications in quantum information
science and technology, review1 ; review2 ; review3 ; review4  quantum
sensing, Jing24  and quantum networks. LiPRXQ

In particular, the coherent couplings among magnons, microwave cavity photons,
and vibration phonons form the system of cavity magnomechanics. CMMreview  In
this system, the magnon mode couples to the cavity (mechanical) mode via the
magnetic dipole (magnetostrictive) interaction. Theoretical studies indicate
that the cavity magnomechanical system offers a promising platform for
preparing macroscopic quantum states, such as entangled states Li2018 ;
Li2019A ; Li2021  and squeezed states Li2019B ; WangHF2021 ; LiNSR  of magnons
and phonons, and ground-state cooling of mechanical motion. Ding20 ; Jing21 ;
Asjad2023  The magnetostriction induced mechanical displacement can further
couple to an optical cavity via radiation pressure, which leads to the
formation of the optomagnomechanical (OMM) system. CMMreview ; Fan2023A  The
OMM system can act as a foundation for building highly hybrid systems, e.g.,
by coupling the optical cavity to an ensemble of atoms Fan2023C , and the
magnon mode to a microwave cavity. Fan2023D  Relevant entanglement protocols
have been put forward, Fan2023A ; Fan2023C ; Fan2023D ; YRC ; DiKe  indicating
that the OMM system has many potential applications in quantum information
processing and quantum technology.

Here, we propose to prepare magnonic squeezed states using the novel OMM
system. We note that several proposals have been offered to prepare magnonic
squeezed states. Specifically, they can be achieved by transferring squeezing
from microwaves to magnons, Li2019B  two-tone driving of the magnon mode,
WangHF2021 ; Qian2024  two-tone driving of a qubit coupled to the magnon mode,
Guo2023  and exploiting the anisotropy of the ferromagnet, Kamra2016 ;
Sharma2021  the nonlinearity of the magnetostriction, LiNSR  and the magnon
Kerr effect due to the magnetocrystalline anisotropy, Asjad2023  etc. In this
work, we adopt a different approach by applying the idea of reservoir
engineering to the OMM system, and consider two scenarios depending on whether
the magnomechanical coupling is linear or dispersive. The strong linear
coupling is typically the case where the mechanical frequency is close to the
magnon frequency and occurs for higher-frequency phonon modes (in the GHz
range) associated with thin YIG films, An20 ; Xu21  while the weak dispersive
coupling arises for lower-frequency phonon modes (in the MHz range) being much
lower than the magnon frequency associated with, e.g., macroscopic YIG
spheres. Zhang2016 ; Li2018 ; Potts2021 ; Shen2022  We first generate strong
squeezing of the mechanical mode by driving the optical cavity with a two-tone
laser field . We then show that in both cases mechanical squeezing can be
efficiently transferred to the magnon mode: Stationary magnon squeezing is
obtained in the linear coupling case; while transient squeezing is achieved in
the dispersive coupling case using a two-step protocol.

The article is organized as follows. In Sec. II, we introduce the OMM system
with the linear magnomechanical coupling and provide the Hamiltonian and
Langevin equations. We then present the results of the stationary magnon
squeezing in Sec. III. In Sec. IV, we describe the OMM system with the
dispersive magnomechanical coupling and introduce the two-step operations for
preparing a transient magnonic squeezed state. We then present the results and
discuss the details related to the experimental realization. Finally, we
summarize our findings in Sec. V.

##  II The OMM system with linear magnomechanical interaction

Figure 1: (a) Schematic diagram of the magnon squeezing protocol using the OMM
system with a linear magnomechanical coupling. The mechanical mode ($b$)
couples to an optical cavity mode ($a$) via the radiation-pressure interaction
and to a magnon mode ($m$) via a linear beam-splitter interaction. (b) Mode
and drive frequencies used in the protocol. Strong mechanical squeezing is
achieved by simultaneously driving the cavity with two laser fields at
frequencies $\omega_{\pm}=\omega_{a}\pm\omega_{b}$, and the squeezing is
further transferred to the magnon mode via the phonon-magnon linear coupling.

The OMM system consists of optical cavity photons, ferrimagnetic magnons, and
vibration phonons induced by magnetostriction, as depicted in Fig. 1(a).
Magnons, quanta of spin wave, can couple to phonons via the magnetostrictive
interaction. A linear and strong magnon-phonon coupling becomes dominant when
the mechanical frequency is high and nearly resonant with the magnon
frequency. CMMreview  This is typically the case when the magnomechanical
sample is a thin YIG film, An20 ; Xu21  and the strong magnon-phonon coupling
leads to magnon polarons. Kikkawa  The magnomechanical displacement further
couples to an optical cavity mode via the radiation pressure or photoelastic
effect, COMreview ; Fan2023A  which is a dispersive type interaction.
Therefore, in the situation of a linear magnon-phonon coupling, the
Hamiltonian of the OMM system reads

$\begin{split}H/\hbar=&\sum_{o=a,b,m}\omega_{o}o^{\dagger}o-g_{0}a^{\dagger}a(b+b^{\dagger})+g_{m}(m^{\dagger}b+mb^{\dagger})\\\
&+[(E_{+}e^{-i\omega_{+}t}+E_{-}e^{-i\omega_{-}t})a^{\dagger}+\mathrm{H.c.}],\end{split}$
(1)

where $o=a$, $b$ and $m$ ($o^{\dagger}$) are the annihilation (creation)
operators of the optical cavity, mechanical, and magnon modes, respectively,
and $\omega_{o}$ are the corresponding resonance frequencies. $g_{m}$ denotes
the coupling strength between the magnon mode and the mechanical mode, which
can be strong; $g_{0}$ is the bare optomechanical coupling strength, which is
typically weak, but the effective optomechanical coupling strength $G$ can be
significantly enhanced by driving the cavity with a strong laser field. The
last term is the driving Hamiltonian corresponding to a two-tone laser field
driving the cavity. Such a two-tone driving strategy has been adopted to
prepare squeezed states of the mechanical motion in cavity optomechanics Tan13
; Kronwald2013 ; Wollman2015 ; Pirkkalainen2015 ; Teufel2015  and cavity
magnomechanics. WangHF2021 ; Qian2024  Specifically, Fig. 1(b) shows the
frequencies of the driving fields at the Stokes
($\omega_{-}=\omega_{a}-\omega_{b}$) and anti-Stokes
($\omega_{+}=\omega_{a}+\omega_{b}$) scattering sidebands of the optical
cavity. Here
$E_{\pm}=\sqrt{\kappa_{a}P_{\pm}/\hbar\omega_{\pm}}e^{i\phi_{\pm}}$ represent
the coupling strengths between the driving fields and the cavity mode, with
$\kappa_{a}$ being the dissipation rate of the cavity, and $P_{\pm}$
($\phi_{\pm}$) the powers (phases) of the driving fields.

Including the dissipation and input noise of each mode, we obtain the
following quantum Langevin equations (QLEs) of the system:

$\begin{split}\dot{a}=&-i\omega_{a}a-\frac{\kappa_{a}}{2}a+ig_{0}a(b+b^{\dagger})\\\
&-i(E_{+}e^{-i\omega_{+}t}+E_{-}e^{-i\omega_{-}t})+\sqrt{\kappa_{a}}a_{in},\\\
\dot{b}=&-i\omega_{b}b-\frac{\gamma_{b}}{2}b+ig_{0}a^{\dagger}a-ig_{m}m+\sqrt{\gamma_{b}}b_{in},\\\
\dot{m}=&-i\omega_{m}m-\frac{\kappa_{m}}{2}m-ig_{m}b+\sqrt{\kappa_{m}}m_{in},\end{split}$
(2)

where $\gamma_{b}$ ($\kappa_{m}$) is the dissipation rate of the mechanical
(magnon) mode. $o_{in}$ ($o=a$, $b$, $m$) are the corresponding input noise
operators for the three modes, which are zero-mean and obey the following
correlation functions: $\langle
o_{in}(t)o_{in}^{\dagger}(t^{\prime})\rangle=[N_{o}(\omega_{o})+1]\delta(t-t^{\prime})$
and $\langle
o_{in}^{\dagger}(t)o_{in}(t^{\prime})\rangle=N_{o}(\omega_{o})\delta(t-t^{\prime})$.
Here, $N_{o}(\omega_{o})=[\mathrm{exp}(\hbar\omega_{o}/k_{B}T)-1]^{-1}$ are
the mean thermal excitation number of each mode at bath temperature $T$, with
$k_{B}$ the Boltzmann constant.

Since the cavity is strongly driven by the two-tone laser field, which
validates the linearization treatment, Qian2024  we can write each mode
operator as the sum of its classical average and quantum fluctuation, i.e.,
$o=\langle o\rangle+\delta o$ ($o=a$, $b$, $m$). Substituting these
decompositions into Eq. (2), we obtain two sets of equations for the classical
averages and the quantum fluctuations, respectively. Since the amplitude of
the cavity field is dominant at the two driving frequencies $\omega_{\pm}$, we
assume that $\langle a\rangle\simeq
a_{+}e^{-i\omega_{+}t}+a_{-}e^{-i\omega_{-}t}$. Kronwald2013 ; Qian2024  In
typical optomechanical experiments, the cavity frequency shift
$2g_{0}\mathrm{Re}\langle b\rangle$ due to the radiation pressure is much
smaller than the mechanical frequency $\omega_{b}$. Therefore, the small
frequency shift term $ig_{0}\langle a\rangle(\langle b\rangle+\langle
b\rangle^{\ast})$ can be safely neglected in solving the set of equations for
the classical averages. As a result, we obtain the mean amplitudes

$a_{\pm}=\frac{-iE_{\pm}}{i(\omega_{\mp}-\omega_{a})+\kappa_{a}/2}=\frac{-iE_{\pm}}{\kappa_{a}/2\mp
i\omega_{b}},$ (3)

associated with the two frequency components of the cavity field.

The set of linearized equations for the quantum fluctuations of the system, in
the interaction picture with respect to
$\sum_{o=a,b,m}\hbar\omega_{o}o^{\dagger}o$, is given by

$\begin{split}\delta\dot{a}=&-\frac{\kappa_{a}}{2}\delta
a+i(G_{-}+G_{+}e^{-2i\omega_{b}t})\delta b\\\
&+i(G_{+}+G_{-}e^{2i\omega_{b}t})\delta
b^{\dagger}+\sqrt{\kappa_{a}}a_{in},\\\
\delta\dot{b}=&-\frac{\gamma_{b}}{2}\delta b-ig_{m}\delta
m+i(G_{-}^{\ast}+G_{+}^{\ast}e^{2i\omega_{b}t})\delta a\\\
&+i(G_{+}+G_{-}e^{2i\omega_{b}t})\delta
a^{\dagger}+\sqrt{\gamma_{b}}b_{in},\\\
\delta\dot{m}=&-\frac{\kappa_{m}}{2}\delta m-ig_{m}\delta
b+\sqrt{\kappa_{m}}m_{in},\end{split}$ (4)

where $G_{\pm}=g_{0}a_{\pm}$ are the effective optomechanical coupling
strengths associated with the two driving fields, which are complex in
general, but can be set real by adjusting the phases $\phi_{\pm}$ of the
driving fields. For simplicity, hereafter we consider real and positive
couplings $G_{\pm}$.

Equation (4) is difficult to solve due to the time-dependent oscillating
terms. However, under the conditions of $\kappa_{a(m)}$, $\gamma_{b}$,
$G_{\pm}\ll\omega_{b}$, Qian2024  we can take the rotating-wave approximation
(RWA) and safely neglect those high-frequency oscillating terms. Consequently,
we obtain the QLEs in the quadrature form, which can be cast in the following
matrix form:

$\dot{u}(t)=Au(t)+n(t),$ (5)

where $u(t)=[\delta X_{a}(t),\delta Y_{a}(t),\delta q(t),\delta p(t),\delta
X_{m}(t),\delta Y_{m}(t)]^{T}$, with $\delta X_{k}=(\delta k+\delta
k^{\dagger})/\sqrt{2}$ and $\delta Y_{k}=i(\delta k^{\dagger}-\delta
k)/\sqrt{2}$ ($k=a$,$m$) being the fluctuations of the amplitude and phase
quadratures of the cavity and magnon modes, respectively, and $\delta
q=(\delta b+\delta b^{\dagger})/\sqrt{2}$ and $\delta p=i(\delta
b^{\dagger}-\delta b)/\sqrt{2}$ being the fluctuations of the mechanical
position and momentum. The drift matrix is given by

$A=\begin{pmatrix}-\frac{\kappa_{a}}{2}&0&0&G_{+}-G_{-}&0&0\\\
0&-\frac{\kappa_{a}}{2}&G_{+}+G_{-}&0&0&0\\\
0&G_{+}-G_{-}&-\frac{\gamma_{b}}{2}&0&0&g_{m}\\\
G_{+}+G_{-}&0&0&-\frac{\gamma_{b}}{2}&-g_{m}&0\\\
0&0&0&g_{m}&-\frac{\kappa_{m}}{2}&0\\\
0&0&-g_{m}&0&0&-\frac{\kappa_{m}}{2}\end{pmatrix},$ (6)

and the vector of input noises
$n(t)=[\sqrt{\kappa_{a}}X_{a}^{in}(t),\sqrt{\kappa_{a}}Y_{a}^{in}(t),$
$\sqrt{\gamma_{b}}q^{in}(t),\sqrt{\gamma_{b}}p^{in}(t),\sqrt{\kappa_{m}}X_{m}^{in}(t),\sqrt{\kappa_{m}}Y_{m}^{in}(t)]^{T}$.

Because of the dynamics of the system is linearized and the input noises are
Gaussian, the state of the system at any given time is Gaussian and thereby
can be characterized by a $6\times 6$ covariance matrix (CM) $V$, with its
entries defined as $V_{ij}=\langle u_{i}(t)u_{j}(t)+u_{j}(t)u_{i}(t)\rangle/2$
($i,j=1,2,...,6$). The steady-state CM $V$ can be directly achieved by solving
the Lyapunov equation Vitali2007

$AV+VA^{T}=-D,$ (7)

where the diffusion matrix
$D=\mathrm{diag}[\kappa_{a}/2,\kappa_{a}/2,\gamma_{b}(N_{b}(\omega_{b})+1/2),\gamma_{b}(N_{b}(\omega_{b})+1/2),\kappa_{m}(N_{m}(\omega_{m})+1/2),\kappa_{m}(N_{m}(\omega_{m})+1/2)]$,
which is defined via $D_{ij}\delta(t-t^{\prime})=\langle
n_{i}(t)n_{j}(t^{\prime})+n_{j}(t^{\prime})n_{i}(t)\rangle/2$. Here we assume
$N_{a}(\omega_{a})\simeq 0$, due to the high frequency of the optical mode.

##  III Stationary quantum squeezing of mechanical and magnon modes

Figure 2: The degree of squeezing of the mechanical displacement (a) and
magnon phase quadrature (b) versus $G_{+}/G_{-}$. The Wigner function of the
squeezed mechanical mode (c) and magnon mode (d), with $G_{+}/G_{-}=0.76$. See
text for the other parameters.

The physical mechanism of generating quantum squeezing under the two-tone
driving field can be explained as follows. Under the RWA, the two driving
fields at the blue and red sidebands ($\omega_{\pm}=\omega_{a}\pm\omega_{b}$)
of the optical cavity can respectively activate the optomechanical parametric
down conversion and beam-splitter interactions, and the effective
optomechanical interaction can thus be denoted by $\delta
a^{\dagger}(G_{+}\delta b^{\dagger}+G_{-}\delta b)+\mathrm{H.c.}$, cf., Eq.
(4). Using the standard squeezing transformation, CMMreview ; Kronwald2013  we
introduce the Bogoliubov mode operator for the mechanical mode, $\delta
B=\mathrm{cosh}(r)\delta b+\mathrm{sinh}(r)\delta b^{\dagger}$, where the
squeezing parameter $r$ is defined as
$r=\frac{1}{2}\mathrm{ln}\frac{G_{-}+G_{+}}{G_{-}-G_{+}}$. The effective
optomechanical interaction is thus rewritten as $\tilde{G}\delta
a^{\dagger}\delta B+\mathrm{H.c.}$, where
$\tilde{G}=\sqrt{G_{-}^{2}-G_{+}^{2}}$ denotes the coupling strength between
the cavity mode and the Bogoliubov mode. The above cavity-Bogoliubov beam-
splitter interaction can lead to the Bogoliubov mode being significantly
cooled to its quantum ground state, which corresponds to the squeezed vacuum
state of the mechanical mode. By further exploiting the strong phonon-magnon
beam-splitter coupling, the squeezing can be efficiently transferred from the
mechanical mode to the magnon mode. Consequently, stationary squeezing of the
mechanical and magnon modes can be generated under the continuous driving of
the two-tone laser field.

The degree of squeezing in the dB unit is defined as

$S=-10\ \mathrm{log}_{10}[\langle\delta Q(t)^{2}\rangle/\langle\delta
Q(t)^{2}\rangle_{\mathrm{vac}}],$ (8)

where $\langle\delta Q(t)^{2}\rangle$ denotes the variance of the quadrature
of one specific mode and $\langle\delta
Q(t)^{2}\rangle_{\mathrm{vac}}=\frac{1}{2}$ corresponds to the vacuum
fluctuation. In Figs. 2(a) and 2(b), we show the stationary squeezing of the
mechanical displacement and magnon phase quadrature, using the following
feasible parameters: COMreview ; Carlos2020B
$\omega_{m}/2\pi=\omega_{b}/2\pi=10\ \mathrm{GHz}$, cavity resonant wavelength
$\lambda=1064\ \mathrm{nm}$, $\kappa_{a}/2\pi=10^{2}\ \mathrm{MHz}$,
$\kappa_{m}/2\pi=1\ \mathrm{MHz}$, $\gamma_{b}/2\pi=10^{4}\ \mathrm{Hz}$,
$g_{0}/2\pi=10^{2}\ \mathrm{kHz}$, $g_{m}/2\pi=10\ \mathrm{MHz}$, $T=10\
\mathrm{mK}$, and a fixed $G_{-}/2\pi=15\ \mathrm{MHz}$, corresponding to a
pump power of $P_{-}\simeq 26.41\ \mathrm{mW}$ for the driving field at
frequency $\omega_{-}=\omega_{a}-\omega_{b}$. For the stability reason, we
consider $G_{+}<G_{-}$. Qian2024  Clearly, there exists an optimal ratio of
$G_{+}/G_{-}$ for the maximal mechanical (magnon) squeezing, as a result of
the trade-off between the following two effects: On the one hand, a
sufficiently large difference between the two couplings $G_{+}$ and $G_{-}$ is
required to have a sufficiently large cooling rate $\tilde{G}$ of the
Bogoliubov mode; On the other hand, almost equal couplings $G_{+}\approx
G_{-}$ is needed to obtain an extremely large squeezing of the mechanical
mode. Figures 2(c) and 2(d) illustrate the corresponding Wigner function of
the squeezed mechanical and magnon modes, with an optimal ratio of
$G_{+}/G_{-}$.

Figure 3: The degree of stationary magnon squeezing $S_{m}$ versus
$G_{+}/G_{-}$, at different bath temperatures of $T=0.01$ K (solid), $T=0.5$ K
(dashed), and $T=1.0$ K (dot-dahed). The other parameters are the same as in
Fig. 2(b).

In Fig. 3, we study the degree of squeezing of the magnon mode $S_{m}$ versus
the ratio $G_{+}/G_{-}$ at different bath temperatures. As the temperature
rises, the system suffers more thermal noises and a larger cooling rate
$\tilde{G}$ is thus needed to efficiently cool the Bogoliubov mode to its
ground state, which implies a larger difference between $G_{+}$ and $G_{-}$,
i.e., a decreasing optimal ratio $G_{+}/G_{-}$ (for a fixed $G_{-}$). This, in
turn, reduces the degree of squeezing that can be achieved in the steady
state. Nevertheless, we see that stationary magnon squeezing can be achieved
for bath temperature above 1 K, manifesting the robustness of the protocol
against thermal noises.

##  IV Magnon squeezing with dispersive magnomechanical interaction

The above two sections study stationary magnonic squeezed states based on the
OMM system with a linear magnomechanical coupling. However, for lower-
frequency mechanical modes, e.g., of a sub-millimeter sized YIG sphere
Zhang2016 ; Li2018 ; Potts2021 ; Shen2022 ; ShenZ2022  or a micron sized YIG
bridge, Heyroth2019 ; Fan2023A ; Kansanen2021  the magnomechanical interaction
is dominated by the nonlinear dispersive coupling. CMMreview  In this case,
the OMM system has both nonlinear opto- and magnomechanical interactions. The
latter weak dispersive coupling hinders an efficient phonon-to-magnon
squeezing transfer (which prefers a linear strong coupling as discussed in
Secs. II and III). Nevertheless, we show in this section that transient
moderate squeezing of the magnon mode can still be achieved by using a two-
step protocol, as depicted in Fig. 4(a).

Figure 4: (a) Schematic diagram of the two-step magnon squeezing protocol
using the OMM system with a dispersive magnomechanical coupling. (b) Mode and
drive frequencies used in the two-step protocol. In the first step, the
optical cavity is simultaneously driven by two laser fields at frequencies
$\omega_{\pm}=\omega_{a}\pm\omega_{b}$ to generate stationary squeezing of the
mechanical mode. In the second step, the squeezing is transferred from the
mechanical mode to the magnon mode via the magnon-phonon state-swap
interaction, realized by driving the magnon mode with a red-detuned microwave
field.

The Hamiltonian of the OMM system with a dispersive magnon-phonon coupling
reads

$\begin{split}H/\hbar=&\sum_{o=a,b,m}\omega_{o}o^{\dagger}o-g_{0}a^{\dagger}a(b+b^{\dagger})\\\
&+g^{\prime}_{m}m^{\dagger}m(b+b^{\dagger})+H_{\mathrm{dri}}/\hbar,\end{split}$
(9)

where $g^{\prime}_{m}$ is the bare magnomechanical coupling strength, and the
effective coupling can be greatly enhanced by driving the magnon mode with a
strong microwave field. Li2018  The last term is the driving Hamiltonian,
which is different in the two-step operations.

In the first step, we aim to achieve stationary squeezing of the mechanical
mode by driving the optical cavity with a two-tone laser field at frequencies
$\omega_{\pm}=\omega_{a}\pm\omega_{b}$, as shown in Fig. 4(b). Due to the weak
magnomechanical dispersive coupling and the absence of a microwave drive field
to enhance the coupling, the magnon-phonon coupling is very weak. The magnon
mode is essentially decoupled from the driven optomechanical system.
Therefore, the system under consideration in the first step can be reduced to
a two-mode model without the magnon mode. The Hamiltonian then becomes

$\begin{split}H_{\mathrm{1st}}/\hbar=&\
\omega_{a}a^{\dagger}a+\omega_{b}b^{\dagger}b-g_{0}a^{\dagger}a(b+b^{\dagger})\\\
&+[(E_{+}e^{-i\omega_{+}t}+E_{-}e^{-i\omega_{-}t})a^{\dagger}+\mathrm{H.c.}],\end{split}$
(10)

which gives rise to the following QLEs

$\begin{split}\dot{a}=&-i\omega_{a}a-\frac{\kappa_{a}}{2}a+ig_{0}a(b+b^{\dagger})\\\
&-i(E_{+}e^{-i\omega_{+}t}+E_{-}e^{-i\omega_{-}t})+\sqrt{\kappa_{a}}a_{in},\\\
\dot{b}=&-i\omega_{b}b-\frac{\gamma_{b}}{2}b+ig_{0}a^{\dagger}a+\sqrt{\gamma_{b}}b_{in}.\end{split}$
(11)

Following the linearization treatment and taking the RWA as used in Sec. II,
we obtain the steady-state averages for the optical cavity $\langle
a\rangle\simeq a_{+}e^{-i\omega_{+}t}+a_{-}e^{-i\omega_{-}t}$ and the
mechanical mode $\langle b\rangle\simeq
b_{+}e^{-2i\omega_{b}t}+b_{0}+b_{-}e^{2i\omega_{b}t}$, where

$\begin{split}a_{\pm}=&\ \frac{-iE_{\pm}}{\kappa_{a}/2\mp i\omega_{b}},\
b_{+}=\ \frac{ig_{0}a_{-}^{\ast}a_{+}}{\gamma_{b}/2-i\omega_{b}},\\\ b_{0}=&\
\frac{ig_{0}(|a_{-}|^{2}+|a_{+}|^{2})}{\gamma_{b}/2+i\omega_{b}},\ b_{-}=\
\frac{ig_{0}a_{+}^{\ast}a_{-}}{\gamma_{b}/2+3i\omega_{b}}.\end{split}$ (12)

Similarly, the steady-state CM $V_{\mathrm{ab}}$ in terms of quantum
fluctuations of the optical and mechanical quadratures can be achieved by
solving the Lyapunov equation,
$A^{\prime}V_{\mathrm{ab}}+V_{\mathrm{ab}}{A^{\prime}}^{T}=-D^{\prime}$, where
the drift matrix, in the interaction picture with respect to
$\hbar\omega_{a}a^{\dagger}a+\hbar\omega_{b}b^{\dagger}b$, is given by

$A^{\prime}=\begin{pmatrix}-\frac{\kappa_{a}}{2}&0&0&G_{+}-G_{-}\\\
0&-\frac{\kappa_{a}}{2}&G_{+}+G_{-}&0\\\
0&G_{+}-G_{-}&-\frac{\gamma_{b}}{2}&0\\\
G_{+}+G_{-}&0&0&-\frac{\gamma_{b}}{2}\end{pmatrix},$ (13)

with $G_{\pm}=g_{0}a_{\pm}$, and the diffusion matrix
$D^{\prime}=\mathrm{diag}[\kappa_{a}/2,\kappa_{a}/2,\gamma_{b}(N_{b}(\omega_{b})+1/2),\gamma_{b}(N_{b}(\omega_{b})+1/2)]$.

Figure 5: (a) The degree of squeezing $S_{b}$ of the mechanical mode versus
$G_{+}/G_{-}$. (b) The Wigner function of the squeezed mechanical mode, with
$G_{+}/G_{-}=0.95$. See text for the other parameters.

In Fig. 5(a), we plot the degree of squeezing $S_{b}$ of the mechanical mode
versus $G_{+}/G_{-}$, using experimentally feasible parameters:
$\omega_{b}/2\pi=10^{2}\ \mathrm{MHz}$, cavity resonant wavelength
$\lambda=1064\ \mathrm{nm}$, $\kappa_{a}/2\pi=2\ \mathrm{MHz}$,
$\gamma_{b}/2\pi=10^{2}\ \mathrm{Hz}$, $g_{0}/2\pi=1\ \mathrm{kHz}$, and
$T=10\ \mathrm{mK}$. We take $G_{-}/2\pi=0.3\ \mathrm{MHz}$ corresponding to a
power $P_{-}\simeq 0.53$ mW of the driving field at frequency
$\omega_{-}=\omega_{a}-\omega_{b}$. A strong squeezing (13 dB) of the
mechanical displacement is achieved at the optimal ratio $G_{+}/G_{-}\simeq
0.95$, corresponding to a power $P_{+}\simeq 0.48$ mW of the driving field at
frequency $\omega_{+}=\omega_{a}+\omega_{b}$. Figure 5(b) shows the
corresponding Wigner function of the squeezed mechanical mode.

When the stationary and strong mechanical squeezing is achieved in the first
step, we then switch off the two driving fields. After a short period, in the
time scale of $\kappa_{a}^{-1}\ll\tau_{0}\ll\gamma_{b}^{-1}$, when the cavity
photons completely dissipate, while the mechanical state remains almost
unchanged, we then apply a red-detuned microwave field to drive the magnon
mode (Fig. 4(b)) to activate the magnomechanical anti-Stokes scattering, which
realizes the magnon-phonon state-swap interaction. Li2019B  This subsequent
state-swap operation can transfer the squeezing from the mechanical mode to
the magnon mode. Note that in the second step, the system can also reduce to a
two-mode model, because the cavity photons already died out before applying
the magnon drive. The optical cavity is essentially decoupled from the
magnomechanical system, and the Hamiltonian of the magnomechanical system is
given by

$\begin{split}H_{\mathrm{2nd}}/\hbar=&\
\omega_{m}m^{\dagger}m+\omega_{b}b^{\dagger}b+g^{\prime}_{m}m^{\dagger}m(b+b^{\dagger})\\\
&+(\Omega m^{\dagger}e^{-i\omega_{0}t}+\mathrm{H.c.}).\end{split}$ (14)

The last term is driving Hamiltonian, where
$\Omega=\frac{\sqrt{5}}{4}\gamma\sqrt{N}B_{d}e^{i\phi_{0}}$ is the Rabi
frequency, characterizing the coupling strength between the magnon mode and
the driving magnetic field with amplitude $B_{d}$, frequency $\omega_{0}$, and
phase $\phi_{0}$. Li2018  Here $\gamma$ is the gyromagnetic ratio and $N$ is
the total number of spins in the YIG crystal.

Since the two-tone laser drive is switched off (the mechanism of generating
persistent squeezing is absent), the mechanical squeezing gradually reduces
due to its interaction with the thermal bath. Therefore, in the second step
the magnon mode will not be prepared in a stationary squeezed state, but we
shall implement a dynamical squeezing transfer from mechanics to magnons. This
is realized by driving the magnon mode with a microwave field with limited
duration, and the drive is switched off when the magnon mode evolves into a
maximally squeezed state. Under the driving field, the QLEs in the interaction
picture with respect to $\hbar\omega_{0}m^{\dagger}m$ are given by

$\begin{split}\dot{m}=&-i\Delta_{m}m-\frac{\kappa_{m}}{2}m-ig^{\prime}_{m}m(b+b^{\dagger})-i\Omega+\sqrt{\kappa_{m}}m_{in},\\\
\dot{b}=&-i\omega_{b}b-\frac{\gamma_{b}}{2}b-ig^{\prime}_{m}m^{\dagger}m+\sqrt{\gamma_{b}}b_{in},\end{split}$
(15)

where the detuning $\Delta_{m}=\omega_{m}-\omega_{0}$ is set to be
$\Delta_{m}=\omega_{b}$ (Fig. 4(b)) to activate the magnomechanical state-swap
interaction.

Similarly, the above QLEs can be separated into two sets of equations, which
describe the time evolution of the averages and fluctuations, respectively.
Specifically, the dynamical solutions of the averages $\langle m\rangle_{t}$
and $\langle b\rangle_{t}$ can be obtained by numerically solving the
equations

$\begin{split}\dot{\langle m\rangle}_{t}=&-i\Delta_{m}\langle
m\rangle_{t}-\frac{\kappa_{m}}{2}\langle m\rangle_{t}-ig_{m}\langle
m\rangle_{t}(\langle b\rangle_{t}+\langle b\rangle_{t}^{\ast})-i\Omega,\\\
\dot{\langle b\rangle}_{t}=&-i\omega_{b}\langle
b\rangle_{t}-\frac{\gamma_{b}}{2}\langle b\rangle_{t}-ig_{m}\langle
m\rangle_{t}^{\ast}\langle m\rangle_{t},\end{split}$ (16)

with initial conditions $\langle m\rangle_{t=0}$ and $\langle b\rangle_{t=0}$.
The dynamical solutions of the averages are necessary because they will affect
the dynamics of the quantum fluctuations, e.g., the effective coupling
strength $G_{m_{t}}=g^{\prime}_{m}\langle m\rangle_{t}$ in the linearized
equations for quantum fluctuations is determined by $\langle m\rangle_{t}$.

Moving to another interaction picture by introducing the slowly moving
operators with tildes: $\delta\tilde{m}=\delta me^{i\tilde{\Delta}_{m}t}$ and
$\delta\tilde{b}=\delta be^{i\omega_{b}t}$, Li2021  the QLEs for the quantum
fluctuations are given by (with the tildes being dropped for simplicity)

$\begin{split}\delta\dot{m}=&-\frac{\kappa_{m}}{2}\delta m-iG_{m_{t}}(\delta
b+\delta b^{\dagger}e^{2i\omega_{b}t})+\sqrt{\kappa_{m}}m_{in},\\\
\delta\dot{b}=&-\frac{\gamma_{b}}{2}\delta b-i(G_{m_{t}}^{\ast}\delta
m+G_{m_{t}}\delta
m^{\dagger}e^{2i\omega_{b}t})+\sqrt{\gamma_{b}}b_{in},\end{split}$ (17)

where $\tilde{\Delta}_{m}=\Delta_{m}+g^{\prime}_{m}(\langle
b\rangle_{t}+\langle b\rangle_{t}^{\ast})\simeq\omega_{b}$ is the effective
detuning and $G_{m_{t}}=g^{\prime}_{m}\langle m\rangle_{t}$ is the time-
dependent effective coupling strength. Under the conditions of $\kappa_{m}$,
$\gamma_{b}$, $G_{m_{t}}\ll\omega_{b}$, we can take the RWA and neglect the
nonresonant oscillating terms in Eq. (17), and consequently obtain

$\begin{split}\delta\dot{m}=&-\frac{\kappa_{m}}{2}\delta m-iG_{m_{t}}\delta
b+\sqrt{\kappa_{m}}m_{in},\\\ \delta\dot{b}=&-\frac{\gamma_{b}}{2}\delta
b-iG_{m_{t}}^{\ast}\delta m+\sqrt{\gamma_{b}}b_{in},\end{split}$ (18)

which clearly reveal a magnomechanical beam-splitter interaction, and it
enables the transfer of the squeezing from mechanics to magnons. The above
equations can be written in the quadrature form, and cast in the following
matrix form:

$\dot{u}_{\mathrm{mb}}(t)=A^{\prime\prime}(t)u_{\mathrm{mb}}(t)+n^{\prime\prime}(t),$
(19)

where $u_{\mathrm{mb}}(t)=[\delta X_{m}(t),\delta Y_{m}(t),\delta q(t),\delta
p(t)]^{T}$ is the vector of quadrature fluctuations,
$n^{\prime\prime}(t)=[\sqrt{\kappa_{m}}X_{m}^{in}(t),\sqrt{\kappa_{m}}Y_{m}^{in}(t),$
$\sqrt{\gamma_{b}}q^{in}(t),\sqrt{\gamma_{b}}p^{in}(t)]^{T}$ is the vector of
input noises, and

$A^{\prime\prime}(t)=\begin{pmatrix}-\frac{\kappa_{m}}{2}&0&\mathrm{Im}G_{m_{t}}&\mathrm{Re}G_{m_{t}}\\\
0&-\frac{\kappa_{m}}{2}&-\mathrm{Re}G_{m_{t}}&\mathrm{Im}G_{m_{t}}\\\
-\mathrm{Im}G_{m_{t}}&\mathrm{Re}G_{m_{t}}&-\frac{\gamma_{b}}{2}&0\\\
-\mathrm{Re}G_{m_{t}}&-\mathrm{Im}G_{m_{t}}&0&-\frac{\gamma_{b}}{2}\end{pmatrix}$
(20)

is the time-dependent drift matrix. The squeezing of the magnon mode can be
evaluated by solving the dynamical CM $V_{\mathrm{mb}}(t)$ of the
magnomechanical system, via numerically solving

$\dot{V}_{\mathrm{mb}}(t)=A^{\prime\prime}(t)V_{\mathrm{mb}}(t)+V_{\mathrm{mb}}(t)A^{\prime\prime}(t)^{T}+D^{\prime\prime},$
(21)

with the initial CM $V_{\mathrm{mb}}(0)$ at the beginning of the second step
(the time the magnon drive is applied), where
$D^{\prime\prime}=\mathrm{diag}[\kappa_{m}(N_{m}(\omega_{m})+1/2),\kappa_{m}(N_{m}(\omega_{m})+1/2),\gamma_{b}(N_{b}(\omega_{b})+1/2),\gamma_{b}(N_{b}(\omega_{b})+1/2)]$
is the diffusion matrix. To check the validity of the RWA, we also calculate
$V_{\mathrm{mb}}(t)$ by numerically solving Eq. (17) without taking the RWA
for comparison.

Figure 6: Time evolution of the degree of magnon squeezing $S_{m}$ in the
second step for different switch-off times of the two-tone drive: (a)
$\omega_{b}t^{\prime}=0.5\pi$, (b) $\omega_{b}t^{\prime}=0$, (c)
$\omega_{b}t^{\prime}=0.2\pi$, and (d) $\omega_{b}t^{\prime}=0.92\pi$, with
the RWA (blue curves) and without the RWA (pink curves). See text for the
other parameters.

Figure 6 shows the time evolution of the variance of the magnon phase
quadrature $\langle\delta Y_{m}(t)^{2}\rangle$ starting from the time when the
magnon drive is applied. We adopt experimentally feasible parameters:
$\omega_{m}/2\pi=10$ GHz, $\omega_{b}/2\pi=10^{2}$ MHz, $\kappa_{m}/2\pi=1$
MHz, $\gamma_{b}/2\pi=10^{2}$ Hz, $g^{\prime}_{m}=10$ Hz,
$\Delta_{m}=\omega_{b}$, $T=10$ mK, and $P\simeq 4.93$ mW for a $10\times
3.0\times 1.5\ \mu m^{3}$ YIG microbridge. Fan2023A  Figure 6 is obtained by
numerically solving Eq. (16) with the initial conditions $\langle
m\rangle_{t=0}=0$ and $\langle b\rangle_{t=0}$ and solving Eq. (21) with the
initial $V_{\mathrm{mb}}(0)=[V_{m}^{\mathrm{th}},0;\ 0,\
V_{b}^{\mathrm{opt}}]$. Here,
$V_{m}^{\mathrm{th}}=\mathrm{diag}[N_{m}(\omega_{m})+\frac{1}{2},N_{m}(\omega_{m})+\frac{1}{2}]$
denotes that the magnon mode is initially in a thermal state with the mean
thermal excitation number $N_{m}$, and $V_{b}^{\mathrm{opt}}$ is the CM of the
squeezed mechanical mode prepared in the first step (cf., Fig. 5(b)), which
remains almost unchanged in the free evolution process with duration of
$\tau_{0}\ll\gamma_{b}^{-1}$. Specifically, we consider
$\kappa_{a}^{-1}\ll\tau_{0}=\frac{4\pi}{\kappa_{a}}\ll\gamma_{b}^{-1}$, which
is long enough for cavity photons to completely decay.

We note that the initial condition $\langle b\rangle_{t=0}$ used for solving
Eq. (16) depends on the time the two-tone drive is switched off. This can be
seen from the results of Eq. (12), where the steady-state averages of the
cavity and mechanical modes are periodically oscillating, a typical feature
with a two-tone drive. Therefore, different turn-off times (of the two-tone
drive) give different initial conditions for solving the equations governing
the free evolution process with duration $\tau_{0}$, i.e., the QLEs (11)
without the driving terms. This further leads to different initial conditions
$\langle b\rangle_{t=0}$ for solving Eq. (16) in the second step. To study the
influence of the switch-off time, we show in Fig. 6 the magnon squeezing for
four different turn-off times within a mechanical period
$\omega_{b}t^{\prime}=0$, $0.2\pi$, $0.5\pi$ and $0.92\pi$ after the
mechanical mode reaches the steady state under the two-tone drive. Clearly,
the dynamical evolution of the magnon squeezing is almost unaffected. To
verify the validation of the RWA used for plotting Fig. 6 (blue curves), we
also compare the results with those obtained without taking the RWA (pink
curves), i.e., by solving Eq. (17). Obviously, the RWA is a very good
approximation under the parameters of Fig. 6.

Figure 6 indicates that a dynamical magnon mode squeezing of $S_{m}^{\rm
max}=2.3$ dB can be achieved in the second step. The microwave driving field
can be turned off when the squeezing reaches its maximum, and the transient
magnon squeezing can be further transferred to a microwave field using their
state-swap interaction Li2019B , which is useful for quantum information
processing or quantum metrology tasks.

##  V Conclusions

We have presented two protocols for preparing magnonic squeezed states in two
configurations of an OMM system by exploiting reservoir engineering. In the
linear (dispersive) magnomechanical coupling case, stationary (transient)
magnonic squeezed states can be generated. The magnon squeezing is transferred
from the squeezed mechanical mode that is produced via two-tone driving of the
optical cavity. The phonon-magnon state-swap interaction is responsible for
the squeezing transfer, which is an intrinsic interaction in the linear
coupling case and can be effectively activated in the dispersive coupling case
by applying a red-detuned magnon drive. The magnonic squeezed state of a
large-size YIG sphere is a macroscopic quantum state and thus useful in
testing the limits of quantum mechanics. The magnonic squeezed states also
find broad applications in quantum information processing, quantum metrology,
and quantum sensing based on magnons.

## Acknowledgements

We thank J. Cheng and H. Qian for useful discussions. This work has been
supported by National Key Research and Development Program of China (Grant No.
2022YFA1405200) and National Natural Science Foundation of China (Grant No.
92265202).

## References

  * (1) D. Lachance-Quirion, Y. Tabuchi, A. Gloppe, K. Usami, and Y. Nakamura, "Hybrid quantum systems based on magnonics," Appl. Phys. Express 12, 070101 (2019).
  * (2) Y. Li, W. Zhang, V. Tyberkevych, W.-K. Kwok, A. Hoffmann, and V. Novosad, "Hybrid magnonics: Physics, circuits, and applications for coherent information processing," J. Appl. Phys. 128, 130902 (2020).
  * (3) H. Y. Yuan, Y. Cao, A. Kamra, R. A. Duine, and P. Yan, "Quantum magnonics: When magnon spintronics meets quantum information science," Phys. Rep. 965, 1-74 (2022).
  * (4) B. Z. Rameshti, S. V. Kusminskiy, J. A. Haigh, K. Usami, D. Lachance-Quirion, Y. Nakamura, C.-M. Hu, H. X. Tang, G. E. W. Bauer, and Y. M. Blanter, "Cavity magnonics," Phys. Rep. 979, 1-61 (2022).
  * (5) H. Huebl, C. W. Zollitsch, J. Lotze, F. Hocke, M. Greifenstein, A. Marx, R. Gross, and S. T. B. Goennenwein, "High Cooperativity in Coupled Microwave Resonator Ferrimagnetic Insulator Hybrids," Phys. Rev. Lett. 111, 127003 (2013).
  * (6) Y. Tabuchi, S. Ishino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, "Hybridizing Ferromagnetic Magnons and Microwave Photons in the Quantum Limit," Phys. Rev. Lett. 113, 083603 (2014).
  * (7) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, "Strongly Coupled Magnons and Cavity Microwave Photons," Phys. Rev. Lett. 113, 156401 (2014).
  * (8) A. Osada, R. Hisatomi, A. Noguchi, Y. Tabuchi, R. Yamazaki, K. Usami, M. Sadgrove, R. Yalla, M. Nomura, and Y. Nakamura, "Cavity Optomagnonics with Spin-Orbit Coupled Photons," Phys. Rev. Lett. 116, 223601 (2016).
  * (9) X. Zhang, N. Zhu, C.-L. Zou, and H. X. Tang, "Optomagnonic Whispering Gallery Microresonators," Phys. Rev. Lett. 117, 123605 (2016).
  * (10) J. A. Haigh, A. Nunnenkamp, A. J. Ramsay, and A. J. Ferguson, "Triple-Resonant Brillouin Light Scattering in Magneto-Optical Cavities," Phys. Rev. Lett. 117, 133602 (2016).
  * (11) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, "Cavity magnomechanics," Sci. Adv. 2, e1501286 (2016).
  * (12) J. Li, S.-Y. Zhu, and G. S. Agarwal, "Magnon-Photon-Phonon Entanglement in Cavity Magnomechanics," Phys. Rev. Lett. 121, 203601 (2018).
  * (13) C. A. Potts, E. Varga, V. Bittencourt, S. V. Kusminskiy, and J. P. Davis, "Dynamical Backaction Magnomechanics," Phys. Rev. X 11, 031053 (2021).
  * (14) R.-C. Shen, J. Li, Z.-Y. Fan, Y.-P. Wang, and J. Q. You, "Mechanical Bistability in Kerr-modified Cavity Magnomechanics," Phys. Rev. Lett. 129, 123601 (2022).
  * (15) Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, "Coherent coupling between a ferromagnetic magnon and a superconducting qubit," Science 349, 405 (2015).
  * (16) D. Lachance-Quirion, Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki, and Y. Nakamura, "Resolving quanta of collective spin excitations in a millimeter-sized ferromagnet," Sci. Adv. 3, e1603150 (2017).
  * (17) D. Lachance-Quirion, S. P.Wolski, Y. Tabuchi, S. Kono, K. Usami, and Y. Nakamura, "Entanglement-based single-shot detection of a single magnon with a superconducting qubit," Science 367, 425 (2020).
  * (18) D. Xu, X.-K. Gu, H.-K. Li, Y.-C. Weng, Y.-P. Wang, J. Li, H. Wang, S.-Y. Zhu, and J. Q. You, "Quantum Control of a Single Magnon in a Macroscopic Spin System," Phys. Rev. Lett. 130, 193603 (2023).
  * (19) Q. Zhang, J. Wang, T.-X. Lu, F. Nori, and H. Jing, Quantum Weak Force Sensing with Squeezed Magnomechanics, arXiv:2404.00723 (2024).
  * (20) J. Li, Y.-P. Wang, W.-J. Wu, S.-Y. Zhu, and J. Q. You, "Quantum network with magnonic and mechanical nodes," PRX Quantum 2, 040344 (2021).
  * (21) X. Zuo, Z.-Y. Fan, H. Qian, M.-S. Ding, H. Tan, H. Xiong, and J. Li, "Cavity magnomechanics: from classical to quantum," New J. Phys. 26, 031201 (2024).
  * (22) J. Li and S.-Y. Zhu, "Entangling two magnon modes via magnetostrictive interaction," New J. Phys. 21, 085001 (2019).
  * (23) J. Li and S. Gröblacher, "Entangling the vibrational modes of two massive ferromagnetic spheres using cavity magnomechanics", Quantum Sci. Technol. 6, 024005 (2021).
  * (24) J. Li, S.-Y. Zhu, and G. S. Agarwal, "Squeezed states of magnons and phonons in cavity magnomechanics," Phys. Rev. A 99, 021801 (2019).
  * (25) W. Zhang, D.-Y. Wang, C.-H. Bai, T. Wang, S. Zhang, and H.-F. Wang, "Generation and transfer of squeezed states in a cavity magnomechanical system by two-tone microwave fields," Opt. Express 29, 11773 (2021).
  * (26) J. Li, Y.-P. Wang, J.-Q. You, and S.-Y. Zhu, "Squeezing microwaves by magnetostriction," Natl. Sci. Rev. 10, nwac247 (2023).
  * (27) M.-S. Ding, L. Zheng, and C. Li, "Ground-state cooling of a magnomechanical resonator induced by magnetic damping," J. Opt. Soc. Am. B 37, 627-634 (2020).
  * (28) T.-X. Lu, H. Zhang, Q. Zhang, and H. Jing, "Exceptional-point-engineered cavity magnomechanics," Phys. Rev. A 103, 063708 (2021).
  * (29) M. Asjad, J. Li, S.-Y. Zhu, and J. Q. You, "Magnon squeezing enhanced ground-state cooling in cavity magnomechanics," Fundamental Res. 3, 3-7 (2023).
  * (30) Z. Fan, H. Qian, and J. Li, "Stationary optomagnonic entanglement and magnon-to-optics quantum state transfer via opto-magnomechanics," Quantum Sci. Technol. 8, 015014 (2023).
  * (31) Z.-Y. Fan, H. Qian, X. Zuo, and J. Li, "Entangling ferrimagnetic magnons with an atomic ensemble via optomagnomechanics," Phys. Rev. A 108, 023501 (2023).
  * (32) Z.-Y. Fan, L. Qiu, S. Gröblacher, and J. Li, "Microwave-Optics Entanglement Via Cavity Optomagnomechanics," Laser Photonics Rev. 17, 2200866 (2023).
  * (33) Y.-X. Luo, L.-J. Cong, Z.-G. Zheng, H.-Y. Liu, Y. Ming, and R.-C. Yang, "Entanglement enhancement and EPR steering based on a PT-symmetric-like cavity-opto-magnomechanical hybrid system," Opt. Express 31, 34764 (2023)
  * (34) K. Di, X. Wang, H. Xia, Y. Zhao, Y. Liu, A. Cheng, and J. Du, "Macroscopic entanglement between ferrimagnetic magnons and atoms via crossed optical cavities," Opt. Lett. 49, 2878 (2024).
  * (35) H. Qian, X. Zuo, Z.-Y. Fan, J. Cheng, and J. Li, "Strong squeezing of microwave output fields via reservoir-engineered cavity magnomechanics", Phys. Rev. A 109, 013704 (2024).
  * (36) Q. Guo, D. Xu, J. Cheng, H. Tan, and J. Li, "Magnon squeezing by two-tone driving of a qubit in cavity-magnon-qubit systems," Phys. Rev. A 108, 063703 (2023).
  * (37) A. Kamra and W. Belzig, "Super-Poissonian shot noise of squeezed-magnon mediated spin transport," Phys. Rev. Lett. 116, 146601 (2016).
  * (38) S. Sharma, V. A. S. V. Bittencourt, A. D. Karenowska, and S. V. Kusminskiy, "Spin cat states in ferromagnetic insulators," Phys. Rev. B 103, L100403 (2021).
  * (39) K. An, A. N. Litvinenko, R. Kohno, A. A. Fuad, V. V. Naletov, L. Vila, U. Ebels, G. de Loubens, H. Hurdequint, N. Beaulieu, J. Ben Youssef, N. Vukadinovic, G. E. W. Bauer, A. N. Slavin, V. S. Tiberkevich, and O. Klein, "Coherent long-range transfer of angular momentum between magnon Kittel modes by phonons," Phys. Rev. B 101, 060407 (2020).
  * (40) J. Xu, C. Zhong, X. Zhou, X. Han, D. Jin, S. K. Gray, L. Jiang, and X. Zhang, "Coherent Pulse Echo in Hybrid Magnonics with Multimode Phonons," Phys. Rev. Appl. 16, 024009 (2021)
  * (41) T. Kikkawa, K. Shen, B. Flebus, R. A. Duine, K. Uchida, Z. Qiu, G. E. W. Bauer, and E. Saitoh, "Magnon Polarons in the Spin Seebeck Effect," Phys. Rev. Lett. 117, 207203 (2016).
  * (42) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, "Cavity optomechanics," Rev. Mod. Phys. 86, 1391 (2014).
  * (43) H. Tan, G. Li, and P. Meystre, "Dissipation-driven two-mode mechanical squeezed states in optomechanical systems," Phys. Rev. A 87, 033829 (2013).
  * (44) A. Kronwald, F. Marquardt, and A. A. Clerk, "Arbitrarily large steady-state bosonic squeezing via dissipation," Phys. Rev. A 88, 063833 (2013).
  * (45) E. E. Wollman, C. U. Lei, A. J. Weinstein, J. Suh, A. Kronwld, F. Marquardt, A. A. Clerk, and K. C. Schwab, "Quantum squeezing of motion in a mechanical resonator," Science 349, 952 (2015).
  * (46) F. Lecocq, J. B. Clark, R. W. Simmonds, J. Aumentado, and J. D. Teufel, "Quantum Nondemolition Measurement of a Nonclassical State of a Massive Object," Phys. Rev. X 5, 041037 (2015).
  * (47) J.-M. Pirkkalainen, E. Damskägg, M. Brandt, F. Massel, and M. A. Sillanpää, "Squeezing of Quantum Noise of Motion in a Micromechanical Resonator," Phys. Rev. Lett. 115, 243601 (2015).
  * (48) D. Vitali, S. Gigan, A. Ferreira, H. R. Böhm, P. Tombesi, A. Guerreiro, V. Vedral, A. Zeilinger, and M. Aspelmeyer, "Optomechanical Entanglement between a Movable Mirror and a Cavity Field," Phys. Rev. Lett. 98, 030405 (2007).
  * (49) C. Gonzalez-Ballestero, D. Hümmer, J. Gieseler, and O. Romero-Isart, "Theory of quantum acoustomagnonics and acoustomechanics with a micromagnet", Phys. Rev. B 101, 125404 (2020).
  * (50) Z. Shen, G. Xu, M. Zhang, Y. Zhang, Y.Wang, C. Chai, C. Zou, G. Guo, and C. Dong, "Coherent Coupling between Phonons, Magnons, and Photons," Phys. Rev. Lett. 129, 243601 (2022).
  * (51) F. Heyroth, C. Hauser, P. Trempler, P. Geyer, F. Syrowatka, R. Dreyer, S. G. Ebbinghaus, G. Woltersdorf, and G. Schmidt, "Monocrystalline Freestanding Three-Dimensional Yttrium-Iron-Garnet Magnon Nanoresonators," Phys. Rev. Appl. 12, 054031 (2019).
  * (52) K. S. U. Kansanen, C. Tassi, H. Mishra, M. A. Sillanpää, and T. T. Heikkilä, "Magnomechanics in suspended magnetic beams," Phys. Rev. B 104, 214416 (2021).


# Perception-Inspired Graph Convolution for Music Understanding Tasks

Emmanouil Karystinaios1 Francesco Foscarin1,2 Gerhard Widmer1,2 1Institute of
Computational Perception, Johannes Kepler University Linz, Austria
2LIT AI Lab, Linz Institute of Technology, Austria
{emmanouil.karystinaios, francesco.foscarin<EMAIL_ADDRESS>
###### Abstract

We propose a new graph convolutional block, called MusGConv, specifically
designed for the efficient processing of musical score data and motivated by
general perceptual principles. It focuses on two fundamental dimensions of
music, pitch and rhythm, and considers both relative and absolute
representations of these components. We evaluate our approach on four
different musical understanding problems: monophonic voice separation,
harmonic analysis, cadence detection, and composer identification which, in
abstract terms, translate to different graph learning problems, namely, node
classification, link prediction, and graph classification. Our experiments
demonstrate that MusGConv improves the performance on three of the
aforementioned tasks while being conceptually very simple and efficient. We
interpret this as evidence that it is beneficial to include perception-
informed processing of fundamental musical concepts when developing graph
network applications on musical score data. All code and models are released
on https://github.com/manoskary/musgconv.

##  1 Introduction

Music data can be represented in computer applications in multiple formats,
two popular ones being audio and symbolic representations. The first encodes a
measure of the pressure/intensity of the sound wave over time, while the
second explicitly encodes discrete musical events such as notes and rests (see
Foscarin et al. (2022) for an overview of different symbolic music formats).
Due to this higher-level information, symbolic representations are generally
considered better inputs/outputs for music analysis and generation tasks in
the Music Information Research (MIR) field. Moreover, any musical task that
starts from a musical score (sheet music) or from MIDI files is naturally in
the symbolic domain.

In the MIR literature, symbolically encoded music is typically handled with
techniques heavily inspired by computer vision (CV) or natural language
processing (NLP) research. In the first case, music is represented in a so-
called ”piano roll” format (first developed in MIDI sequencers) and treated as
a raster image, with the X axis being time, and the Y axis pitch. The pitch
values are typically encoded as MIDI pitch, i.e., with integers in $[0-127]$,
which covers all notes that can be played by common well-tempered instruments
(a range that is larger than the range of the piano, with 88 notes), and the
time resolution is a parameter that is usually set to the expected shortest
note duration. In its simpler version, each element of the 2D matrix is set to
1 if a note is sounding at the corresponding pitch and time, or 0 otherwise.
The downside of this approach is that it creates a very large and very sparse
input matrix since only a few notes will play at any time. The other common
approach is treating music with sequential models from NLP research. Although
different tokenization techniques have been proposed Fradet et al. (2021), it
is easy to argue that music does not fit well into strictly sequential models,
since more than one note can sound simultaneously, notes can partially
overlap, and generally, the pitch-temporal relations between notes hold
important musical information.

A few recent works Hsiao and Su (2021); Jeong et al. (2019); Karystinaios and
Widmer (2022); Karystinaios et al. (2023); Zhang et al. (2023) started to
explore the use of graphs, and Graph Neural Networks (GNNs), to represent and
process symbolic music. Significant advances have been reported on various
musical tasks Karystinaios et al. (2023); Karystinaios and Widmer (2023).
However, the components that these papers use are “borrowed” from GNN research
with other types of data.

We argue that this can lead to suboptimal results. As a solution, we design
MusGConv, a new graph convolutional block specifically dedicated to music
data, that is founded on fundamental music perception principles. We test our
approach on four musical tasks: voice separation, composer classification,
Roman numeral analysis, and cadence detection. This selection of tasks allows
us to cover three major graph neural network classes of problems: graph
classification, link prediction, and node classification. We compare our
results with the state-of-the-art graph models for these tasks and show that
the use of MusGConv leads to better results overall. Moreover, its simple
design enables this performance boost without adding any additional
computational cost.

##  2 Perceptual and Modeling Considerations

We base our research on the perceptual principles of two fundamental musical
dimensions: pitch and rhythm. Cognitive studies show that people are not very
sensitive to the absolute pitch of individual notes and perceive mainly the
distance between pitches Deutsch (2013). Thus, we perceive the same musical
pattern if it is shifted higher or lower in frequency. This is called relative
pitch perception and has been formalized in music theory through measures of
pitch distance called intervals, which are the basis of all musical concepts
that involve combinations of pitches, such as chords and harmony. The position
(and duration) of notes in time is also not important in itself but only
relative to the position (and duration) of the other notes. This temporal
organization is called rhythm and is typically composed of patterns that tend
to be periodic and organized at different hierarchical levels.

While these principles are simple, producing relative features to use as input
for deep learning systems is not an easy task. For the pitch representation,
considering intervals between consecutive notes instead of absolute pitches,
as proposed for example in the IDyOM framework Pearce (2018, 2005), poses the
problem of defining a note order. This is trivial for monophonic melodies but
becomes problematic for polyphonic music, where multiple notes can fully or
partially overlap in time. Figure 1 (c) exemplifies a possible ordering rule,
but such rules may not be generally valid for different pieces and contexts or
in general yield very different interval sequences for very similar patterns
(e.g., if we remove the F), thus not helping with learning a general
representation. Inserting interval information after the music has been
tokenised in a sequential representation Kermarec et al. (2022), creates
similar ordering problems.

Another strategy is to rely on music theory and consider pitch distances
relative to the fundamental note defined by the key of the piece (see Fig. 1
(b). This has a musicological limitation, as the key signature of a piece may
not change when there are temporary modulations to other keys, and a practical
one, as keys are not always notated in musical datasets and MIDI files. Not to
mention all the music that falls outside the classic tonal framework, for
which the definition of a key is not even meaningful. The most common
alternative, e.g., Nakamura et al. (2016); Foscarin et al. (2021), has been
the use of data augmentation via transposition, assuming that if the network
sees transpositions of the same piece, it will learn patterns that are general
across all transpositions. However, this is far from ideal since the network
will need to store similar patterns redundantly, thus making inefficient use
of network capacity and drastically increasing training time.

For this reason, we believe that designing relative features to input to our
system is not a viable option for general music modelling. We explore a
different path, which is to customize the working mechanism of our network to
take into account the relative perceptual properties of music through a
dedicated message-passing mechanism that computes pairwise pitch and time
representations. This is weakly related to recent work on audio
representations of music that aims at learning transposition-invariant
features Lattner et al. (2018); Elowsson and Friberg (2019); Lattner et al.
(2019) and tempo-invariant rhythmic patterns Di Giorgi et al. (2021). Other
related work targets graphs whose edges encode geometric information Atzmon et
al. (2018); Wang et al. (2019a); Satorras et al. (2021), intending to build
representations that are invariant to operations such as translation. The
ideas from these last approaches (with minor modifications) are also
beneficial for musical tasks, but we show that our music-specific approach
outperforms them.

Figure 1: Three alternative representations of note pitches in a musical
excerpt: (a) absolute representation in terms of MIDI pitch; (b) relative
pitch distance (ignoring the octave) in semitones relative to the fundamental
pitch specified by the key signature (here: C); (c) relative pitch distance in
semitones from the closest preceding note; in case of chords the order is
defined from bottom to top.

Up to this point, we have highlighted the importance of relative pitch and
time representations. However, their absolute values could also be important
depending on the task at hand. For example, instruments would peculiarly
change their timbre as they approach very low or very high notes in their
range, which is recognisable by a listener. In tonal music, the absolute pitch
of notes defines a key signature which could be relevant in the composer
classification task (specific keys can have very specific meanings to
composers, and within a musical tradition). The same goes for absolute time
positions, for example, with patterns happening at the beginning or end of a
piece. Finally, the output of our network may need to be an absolute pitch
(for example, in the Roman numeral analysis task we describe later ), and
therefore we need to retain this information in the network.

This need for considering both absolute and relative pitch and time positions
motivates the design of our new convolutional block, to be described in detail
in Section 4 below.

##  3 Graph Approaches to Musical Tasks

In this section, we describe existing graph modelling approaches to the four
musical tasks we use to evaluate our proposal. They all have a common pipeline
which involves building a graph from a given musical score (see Figure 2) and
using a series of convolutional blocks to produce context-aware hidden
representations for each node. We start by describing the graph-building
procedure and a generic graph convolutional block; we then proceed by
detailing the tasks and the specific network components used to target them.

Figure 2: General architecture of our pipeline. The first part that produces
the hidden node representation is common among all tasks; the last module is
task-specific.

###  3.1 Graph from Musical Scores

A graph is defined as a set of nodes $V$ and a set of edges $E$, where each
edge $(u,v)\in E$ connects the nodes $v,u\in V$. We extend this definition by
considering labelled edges which we model as a triple $(u,r,v)$, where $u,v\in
V$ and $r\in\mathcal{R}$ is a relation type. A graph with edges of multiple
types (we don’t use different types of nodes in this work) is called
heterogeneous. Moreover, we consider an attributed graph, i.e., a graph where
each node is described by a set of features, which we group in the columns of
the matrix $X$. Our attributed heterogeneous graph is defined as $G=(V,E,X)$.

We create such a directed graph from a musical score following the work of
Karystinaios et al. Karystinaios et al. (2023). Each node $v\in V$ corresponds
to one and only one note in the musical score. $\mathcal{R}$ includes 4 types
of relations: onset, during, follow, and silence, corresponding, respectively,
to two notes starting at the same time, a note starting while the other is
sounding, a note starting when the other ends, and a note starting after a
time when no note is sounding. The inverse edges for during, follows, and
silence relations are also created.

The feature matrix $X$ is composed of the following features extracted from
each note of the score: the pitch class, i.e., one of the 12 note names (C,
C#, D, D#, etc.), the octave in $[1,\dots,7]$, the note duration, encoded as a
single float value $d\in[0,1]$ computed as the ratio of the note and bar
durations, passed through a tanh function to limit its value and give more
resolution to shorter notes, as proposed by Karystinaios et al. Karystinaios
et al. (2023). For the task of Roman Numeral Analysis and Cadence Detection,
we add additional specialized features to be consistent with the approach in
the literature we consider in our evaluation Karystinaios and Widmer (2023,
2022).

###  3.2 Graph Convolution Operation

We now present a generic graph convolutional block, to simplify the
description of our music-dedicated approach in the next section. Given an
attributed homogeneous graph $G$, a graph convolution block that updates the
representation of node $u$ for layer $l+1$ can be described as:

$\displaystyle\bm{h}^{(l+1)}_{u}$
$\displaystyle=\psi\left(\bm{h}^{(l)}_{u},\underset{v\in\mathcal{N}(u)}{\textrm{aggregate}}\,(\\{\bm{\eta}_{vu}^{(l)}\\})\right),$
(1) $\displaystyle\bm{\eta}_{vu}^{(l)}$
$\displaystyle=\phi\left(\bm{h}^{(l)}_{v},\bm{h}^{(l)}_{u}\right)$ (2)
$\displaystyle\bm{h}^{(l+1)}_{u}$
$\displaystyle=\sigma\left(\bm{h}^{(l+1)}_{u}\right)$ (3)

where $\textrm{aggregate}(\cdot)$ denotes a differentiable, permutation
invariant aggregation function, e.g., sum, mean, etc.; $\phi$ and $\psi$ are
called edge operation and node operation, respectively, and denote
differentiable learnable functions such as concatenation, sum, or
multiplication, followed by a linear transformation; $\sigma$ denotes a non-
linear function, $\mathcal{N}(u)$ denotes the neighbours of $u$;
$\bm{h}_{u}^{(l)}$ is the hidden representation of node $u$ at layer $l$.

Furthermore, if we want to leverage edge features, Equation 2 becomes:

$\bm{\eta}_{vu}^{(l)}=\phi\left(\bm{h}^{(l)}_{v},\bm{h}^{(l)}_{u},\bm{e}_{vu}\right)$
(4)

where, $\bm{e}_{uv}$ are features of the directed edge connecting node $v$
with node $u$.

When $G$ is heterogeneous, the function $\mathcal{N}(u)$ in Equation 1 is
modified as proposed in Schlichtkrull et al. (2018) to return only the
neighbours nodes which are connected with an edge of type $r$. Equations 1 and
3 are then computed $|\mathcal{R}|$ times and the results
$\mathbf{h_{r}}^{(l)}_{u}$ aggregated in an unique node latent representation
$\bm{h}^{(l)}_{u}$ as:

$\bm{h}^{(l)}_{u}=\underset{r\in\mathcal{R}}{\textrm{aggregate}}\left(\\{\mathbf{h_{r}}^{(l)}_{u}\\}\right)$
(5)

where $\textrm{aggregate}(\cdot)$ denotes a differentiable, permutation
invariant aggregation function such as sum or mean.

We now turn to describe the tasks that we will use for evaluation and the
task-specific part of our graph model.

###  3.3 Monophonic Voice Separation

Voice separation is the task of segmenting a symbolic music piece into an
unknown number of individual monophonic note streams according to musical and
perceptual criteria. Duoane and Pardo Duane and Pardo (2009) framed the
problem as a link prediction task in which two notes are linked if they are
consecutive in the same voice. Karystinaios et al. Karystinaios et al. (2023)
proposed a GNN-based model that reached new state-of-the-art results.

Following their approach, we perform this task by adding a link predictor
module, consisting of an MLP which takes as input the hidden representations
of nodes and for each pair of nodes performs a binary classification between
the “linked” and “not-linked” classes.

The evaluation metric is the binary F1 score, i.e. the F1 score for the
positive class which represents the true links in the ground truth.
Karystinaios et al. also consider the F1 score after a postprocessing phase,
but we don’t use it to keep the number of metrics of reasonable size, and
because, as they report, postprocessing increases the metric in a way which
does not always correlate perfectly with the network performance.

###  3.4 Composer Classification

Composer classification from a symbolic musical score is the task of
identifying the composer of the score from a list of composers. In the graph
problem taxonomy, this falls in the category of global graph classification
tasks. Following the work of Zhang et al. Zhang et al. (2023), we perform this
task by adding a global mean pooling layer to our architecture, which averages
the latent representations built by our GNN blocks, followed by an MLP which
predicts probabilities over composer classes. The composer classification
predictions are evaluated in terms of classification accuracy.

###  3.5 Roman Numeral Analysis

Roman numeral analysis is a branch of analytical musicology whose goal is to
infer the underlying harmony and chord progressions from a musical score. The
result is a set of complex composite labels (the Roman numerals) which
annotate music at the onset level, i.e., for every score position that
corresponds to one or more note onsets. Related work Chen et al. (2018);
Micchi et al. (2020); McLeod and Rohrmeier (2021); López et al. (2021) frames
RN analysis as a multi-task classification problem where every label is broken
down into 5 components (degree, inversion, root, key, and quality) which are
predicted by different classifiers in hard-parameter-sharing setting. Of these
components, degree, inversion, and quality are transposition invariant and the
rest depend on the absolute pitches in the input. A recent approach
Karystinaios and Widmer (2023) considers a graph input and obtains new state-
of-the-art results.

Following this work, we perform this task by adding to our general
architecture an onset edge pooling layer that contracts the latent
representations from the note-wise level to the onset-wise, i.e., it creates a
single vector per unique onset. All vectors are then ordered by time position
and fed into a GRU layer whose output is finally used by the aforementioned
MLP classifiers. In the taxonomy of graph problems, RN analysis falls in
between node classification and subgraph classification because of the effect
of the edge pooling layer.

The evaluation score is the so-called Chord Symbol Recall (CSR), i.e., the
ratio of the total duration of segments where prediction equals annotation vs.
the total duration of annotated segments Harte (2010).

###  3.6 Cadence Detection

Cadence detection is a music analysis task that consists of detecting
cadences, i.e. phrase endings with a strong and specific melodic-harmonic
closure effect, in a musical score. Cadences are important both
musicologically and perceptually and it is known that they relate to
particular voicings and chord progressions; however, their automatic
predictions remain particularly challenging due to the high number of
exceptions and corner cases. A recent graph approach Karystinaios and Widmer
(2022) framed the problem as a multiclass node classification scenario, by
predicting the presence of a cadence and its type for each note.

Following that work, we use a graph autoencoder architecture and the latent
synthetic oversampling technique SMOTE Chawla et al. (2002) to balance the
heavily unbalanced class labels. For the same reason, the reported evaluation
metric is the macro F1 score.

##  4 Our Approach: MusGConv

Similarly to previous works, we use stacked graph convolutions to create a
hidden representation of notes that is then used as input for specific music
tasks (see Figure 2). In our proposed approach, we replace the graph
convolutional blocks in the stack with our novel convolutional block (see
Figure 3).

Figure 3: Visualization of update for node $u$ in our MusGConv block
(considering only one edge type), corresponding to Eqns. 11 and 8.

In terms of the notation introduced in the general description of the graph
convolution in Section 3.2, our convolutional block is characterised by two
core contributions: a way to build the edge features $\bm{e}$, and the choice
of the edge operation $\phi$ in Equation 4.

###  4.1 Edge Features Computation

For each edge between nodes $u,v$, we consider three edge features:
$e^{\textrm{onset}}_{vu},e^{\textrm{dur}}_{vu},e^{\textrm{pitch}}_{vu}$ each
of them encoded as a single scalar corresponding to the distance between
onset, duration, and pitch, respectively.

$\displaystyle e^{\textrm{onset}}_{vu}$ $\displaystyle=|on(u)-on(v)|$
$\displaystyle e^{\textrm{dur}}_{vu}$ $\displaystyle=|dur(u)-dur(v)|$ (6)
$\displaystyle e^{\textrm{pitch}}_{vu}$ $\displaystyle=|pitch(u)-pitch(v)|$

This roughly corresponds to the computation of distance in papers that deal
with geometrical data Satorras et al. (2021) except that in their case it is a
multidimensional space distance, while we compute a set of one-dimensional
distances since it makes no sense to mix duration, pitch and onset
information. To keep these values in a convenient numerical range, we
normalise each feature with $\ell_{2}$ normalisation over all edges in a
batch.

Additionally, we inform the network about the pitch-class interval (PCInt),
i.e., the distance between notes without considering the octave, or the
interval direction. This can be seen as a relative version of the chroma
feature, which is commonly used in MIR tasks related to the harmonic content
of the music. This integer in $[0,\dots,11]$ is passed through an embedding
layer, i.e., a learnable look-up table which maps these integers to points
$\bm{e}_{vu}^{\textrm{PCInt}}$ in a continuous multidimensional space.

For each edge, all the aforementioned edge features are concatenated in a
single vector:

$\bm{e}_{uv}^{(0)}=\textrm{cat}\left(e^{\textrm{onset}}_{vu},e^{\textrm{dur}}_{vu},e^{\textrm{pitch}}_{vu},\bm{e}_{vu}^{\textrm{PCInt}}\right)$
(7)

###  4.2 Edge Operation

Our second contribution, needed to properly leverage the edge feature
information, is a new formulation of the message passing paradigm (see Figure
3 for a graphical representation). The new edge operation defines the edge
operation $\phi$ in Equation 4 as follows:

$\bm{\eta}^{(l)}_{vu}=\textrm{cat}\left(\bm{W}^{(l)}_{2}\bm{h}^{(l)}_{v},g_{\mathbf{\Theta}}^{(l)}\left(\bm{e}^{(l)}_{vu}\right)\right)$
(8)

where $g_{\mathbf{\Theta}}$ denotes a simple two-layer MLP with a Relu
activation and layer-wise normalization.

Note how we concatenate the edge features here, while other approaches that
deal with edge features tend to apply a permutation-invariant operation such
as multiplication or addition Satorras et al. (2021); Hu et al. (2019). In
this way, the transposition- and time-invariant information carried by the
edge features is treated in the same way as the node feature when computing
the pairwise representation $\eta_{uv}$. This is inspired by musicological
considerations: we don’t want to weight/modify the absolute representation
according to the relative representation, but rather to just use it as input,
as is done by cognitively plausible musical models Pearce (2018, 2005). This
is similar to what we discuss in Section 1 but, by using MusGConv, we no
longer have the problem of setting a (debatable) order of the input notes,
since for each node, we will consider the relative features according to every
other node connected to it (see Figure 4).

Figure 4: Relative Pitch features $e^{\textrm{pitch}}_{vu}$ for the
highlighted note $u$.

We consider two variants of our system which differ in the edge features which
are passed to layers after the first. The first variant, named MusGConv uses
the absolute difference of the node hidden embeddings, i.e. $\forall l>0$,

$\bm{e}^{(l)}_{vu}=|\bm{h}_{v}^{(l)}-\bm{h}_{u}^{(l)}|$ (9)

The second variant, named MusGConv(+EF), uses the edge hidden representation
from the previous layer as the edge features of the next layer:

$\mathbf{e}_{vu}^{(l+1)}=g_{\Theta}^{(l)}\left(\mathbf{e}_{vu}^{(l)}\right)$
(10)

###  4.3 Node Operation

To complete the explanation of our convolutional block, we specify the
specific computation for the generic node operation $\psi$ introduced in
Equation 1:

$\bm{h}_{u}^{(l+1)}=\bm{W}_{1}^{(l)}\textrm{cat}\left(\bm{h}^{(l)}_{u},\sum_{v\in\mathcal{N}(u)}\bm{\eta}^{(l)}_{vu}\right)\\\
$ (11)

The aggregation function is a sum, and we choose to concatenate the hidden
representation of $u$ with the messages $\eta_{vu}$ from the other nodes. For
heterogeneous convolutions (i.e. Equation 5), we use $\textrm{mean}(\cdot)$ as
the aggregation function. Such choices are not motivated by musical reasoning,
and experimenting with other operations could be an interesting direction, but
are out of the scope of this paper.

##  5 Data

Musical score graph datasets are different from common datasets in the graph-
related literature regarding the number and size of graphs. Node
classification and link prediction datasets often only consist of a single
huge graph, coupled with a sampling strategy to obtain subgraphs to train and
evaluate the Graph Neural Network Hamilton et al. (2017); Zou et al. (2019).
On the other hand, graph classification datasets often consist of a large
number of small graphs, with less than 50 nodes Wu et al. (2021). Musical
score graphs are neither small nor extremely large and may vary significantly
in size; a Bach Chorale may have $\sim 100$ notes whereas a Beethoven Sonata
might have more than $5000$ (with every note corresponding to a node in the
graph).

Moreover, popular sampling strategies, such as node-wise sampling, subgraph
sampling, random-walk sampling, and spectral sampling, may yield musically
problematic note configurations, for example, by segregating notes that are
played at the same time while grouping in the same subgraph notes that are
very far apart. This is problematic, especially for musically-local tasks such
as voice separation or Roman numeral analysis, where the system cannot be
expected to produce a meaningful result, and could even learn to perform the
task in the wrong way if some notes (and onsets for the roman numeral
analysis) are missing.

In the previous works on graph scores that we are considering, the problem is
avoided by training on mini-batches which consist of single pieces. However,
this is not an efficient solution since it always leaves a big part of memory
unused, thus unnecessarily prolonging the training time and decreasing the
variability in the batch. We describe the new sampling mechanism we use in the
following section, and then move on to detailing the datasets used in our
experiments.

###  5.1 Data Sampling

While our nodes can be ordered by multiple features, the organizational aspect
that is most prominent from a perceptual point of view is time. Indeed, people
would still recognize a music piece if it is segmented over the time axis,
while, for example, considering only pitches in a certain interval could lead
to meaningless results. There is also perceptual evidence that the offset time
of a note is much less salient than the onset time, especially for percussive
instruments (including the piano) whose sound naturally decays over time
Klapuri and Davy (2006). Therefore, when we create our graphs from a musical
score, we set the node order first by the absolute time of onset and then by
pitch.

Once this ordering is set (and having defined an index function $ind(\cdot)$
that given a node returns its index in this ordering) we randomly sample a
subgraph of size $K>1$ from a piece with $N$ notes, with the following
procedure. If $N>K$ we select the nodes $u$ with $ind(v)\leq ind(u)\leq
ind(v)+N$, where $v$ is a random node which satisfies the inequality
$ind(v)<K-N$. If $N\leq K$ we select all nodes in the piece. Note that we can
still have the problem of segregating notes from the same chord, but this can
now only happen at the temporal boundaries of our subgraph, limiting its
impact on the network.

We can then create a batch consisting of $B$ subgraphs of at most size $N$.
Our batching mechanism uses the approach of Hamilton et al. Hamilton et al.
(2017), i.e., all graphs are joined together in a single batched graph that
will contain $B$ disjoint subgraphs.

###  5.2 Datasets

We use four distinct datasets for our four tasks.

####  5.2.1 Voice Separation

The Graph Voice Separation dataset was introduced by Karystinaios et al.
Karystinaios et al. (2023). This dataset contains graph data created from five
collections: 370 Bach Chorales, 48 Preludes and 48 Fugues from the Bach Well-
Tempered Clavier (Books I and II), 15 Bach Inventions, 15 Bach Sinfonias, and
210 movements from Haydn String Quartets. It contains in total $726,246$ nodes
and $3,408,679$ edges from $1,054$ unique score graphs. Karystinaios et al.
only test on single collections to understand the differences in performance
for different composing styles. To have a single general performance
indicator, we introduce a new data split that uses $70$% of the data for
training, $10$% for validation, and $20$% for testing. This split preserves
the percentage of pieces in each collection and is independent of the size of
each score graph.

####  5.2.2 Composer Clasification

For composer classification, we use the scores from the DCML corpora
dataset111https://github.com/DCMLab/dcml_corpora. The dataset includes 10
composers, for a total of 419 scores, from where we build score graphs with
collectively $710,240$ nodes and $3,924,655$ edges. We create a random data
split with $70$% of the data for training, $10$% for validation, and $20$% for
testing, which preserves the percentage of composers in each set.

####  5.2.3 Roman Numeral Analysis

The Roman numeral analysis dataset with data augmentation was introduced by
Lopez et al. López et al. (2021). We use their dataset (with augmentations on
the train set). The created graphs collectively contain $8,968,413$ nodes,
$38,390,729$ edges, and $5,096,853$ unique onset positions from $7,988$ scores
(after data augmentation).

####  5.2.4 Cadence Detection

For the Cadence Detection task, we use four distinct annotated datasets, the
Mozart Piano Sonatas Hentschel et al. (2021), Haydn String Quartets Sears et
al. (2017), Mozart String Quartets Allegraud et al. (2019), and Bach WTC
Fugues Giraud et al. (2015). The created graph collectively contains $300,602$
nodes and $1,392,753$ edges from $153$ scores. We create a random data split
with $70$% of the data for training, $10\%$ for validation, and $20$% for
testing.

##  6 Experiments

Our model for each of the four tasks is built on the respective current state-
of-the-art model presented in Section 3. From here on forward, we refer to the
previous state-of-the-art architecture as baseline model, and to the same
architecture with the convolutional blocks replaced with our new ones, as
MusGConv model. These original models serve as the baselines in the following
experiment. We follow the implementations of the publicly available code with
no major modifications, except for the sampling technique that we highlighted
before. For each task, we consider the main evaluation metric proposed in the
original paper (presented in Section 3).

All experiments for a certain task are run with the best hyperparameter
setting specified in the respective papers; this includes 2 GNN layers, and
the convolutional blocks being SageConv Hamilton et al. (2017) for all tasks
except the voice separation where ResGatedGraphConv Bresson and Laurent (2017)
is used. We use a fixed training, validation, and test split for each task,
and every experiment is run 10 times with different NN initialisations on a
single GPU. We used one GTX 1080 Ti GPU with 11 GB of VRAM.

###  6.1 Main Results

| Voice Separation | Composer clf | RNA | Cadence Detection
---|---|---|---|---
| (Link Prediction) | (Graph Classification) | (Node Classification) | (Node Classification)
Previous SOTA Arch | $0.8111\pm 0.058$ | $0.4288\pm 0.031$ | $0.3221\pm 0.010$ | $0.4065\pm 0.011$
MusGConv | $0.8142\pm 0.035$ | $\mathbf{0.5233}\pm 0.032$ | $0.3126\pm 0.015$ | $0.4126\pm 0.016$
MusGConv(+EF) | $\mathbf{0.8436}\pm 0.032$ | $0.3939\pm 0.018$ | $0.3177\pm 0.010$ | $\mathbf{0.4295}\pm 0.009$
Table 1: Experimental comparison with previous SOTA models. The evaluation
metric varies for the different tasks (see the corresponding subsections in
Section 3). Marked in bold are the best results when they are statistically
significant.

The goal of our main experiment is to quantitatively verify whether the use of
MusGConv can improve the results on the four tasks compared to the respective
baseline, i.e., the architectures used in the corresponding state-of-the-art
approaches. Each experiment is run 10 times with the fixed task-specific data
split (as described in the previous section) and different random seeds. We
report ASO significance Del Barrio et al. (2018) with a confidence level
$\alpha=0.05$ and $\epsilon_{\text{min}}<0.1$.

The results, summarised in Table 1, show that MusGConv(+EF) produces
statistically significant better results for the Voice Separation and Cadence
Detection task. In the Composer Classification task, the best-performing model
is MusGConv, while MusGConv(+EF) yields worse results than the baseline. We
suspect that being this a global graph classification task, the edge features
propagation is harder to train, and the simpler mechanism employed in MusGConv
is a better choice. There is no statistically significant difference in the
performance on the Roman Numeral Analysis task. We found two potential
explanations for this behaviour: first, the RNA model is very complex, with
multiple components which could hide the effect of a modification on the graph
encoder. Moreover, the RNA dataset is augmented with transpositions in all
keys, and therefore having transposition-invariant features may only yield
minimal (if any) advantage. An experiment without augmentation is not
possible, since the output of the RNA model depends on the absolute pitches at
the input, and the not-augmented dataset is very small.

Figure 5: Ablation studies.

Regarding execution time, for each task, we compute the ratio between the
average time of the 10 runs for baseline and the 10 runs with MusGConv and
MusGConv(+EF). Aggregated across all tasks, this ratio has an average of $1\pm
0.05$. Thus, the usage of MusGConv has a minimal impact on the final execution
time.

###  6.2 Ablation Studies

We conduct four ablation studies to explore different model variations in
terms of architecture and selection of edge features. The results are reported
in Figure 5, where we also include the results for MusGConv and MusGConv(+EF)
from our main experiment (see previous section) for comparison. For some
variations we consider, there is not a clear winner, meaning that, while the
usage of relative features as edges is beneficial overall, different versions
of our system can perform better on different musical tasks.

##### No concatenation.

We change the feature aggregation function $\phi$ (Eq. 2) from a concatenation
(Eq. 8) to a multiplication, to mimic the operation used in convolutional
blocks that deal with edge features, e.g., Atzmon et al. (2018); Wang et al.
(2019a); Satorras et al. (2021)). The results show that this degrades the
performance for the composer classification and cadence detection tasks while
improving RN analysis and Voice Separation.

##### No Edge Input.

In this study, we ignore our manually built edge features from Section 4.1 and
use node feature difference (see Eq. 9) as edge features for all blocks
(including the first). This is similar to the edge features employed in
EdgeConv Wang et al. (2019b) (though it has to be noted that our message
passing is different from EdgeConv). This degrades the model performance on
all tasks but RN.

##### No PCInt.

We quantify the effect of the PCInt edge feature as a much more music-specific
edge feature, compared to the “more standard” feature distances. We observe
that the usage of this feature improves all tasks.

##### Signed Features.

We evaluate the use of features obtained by removing the absolute value
operation in Eq. 4.1. Indeed this is a more informative input; for example,
with the absolute value we encode the difference between two note durations,
but we lose the information about which is longer. On the other hand, it
increases the input numerical range, and one could argue that the network
already has access to absolute features, and PCInt and edge type encode
similar information. The results show that using signed features is not
beneficial.

##  7 Conclusion and Future Work

This paper has presented a graph convolution block dedicated to music
understanding tasks. Its working mechanism is inspired by perceptual
considerations and permits the propagation of transposition-invariant and
relative timing information in the message-passing process. More generally,
our work enables an elegant and efficient way to use pairwise note features,
which have been long studied and employed in monophonic music, for polyphonic
music processing. Specifically, our approach can be summarized in two core
contributions: pitch and time pairwise functions as edge features, and a new
way of aggregating this information inside the convolutional block. The design
of this block is kept simple to give us a performance advantage without
increasing computation time. We experimentally verify the validity of our
proposition on four rather diverse musical tasks covering three graph-related
problems: graph classification, node classification, and link prediction.

As future work, it would be interesting to evaluate MusGConv on other kind of
music and to investigate the impact of other pairwise note functions, like the
one proposed in the cognitively plausible music model IDyOM Pearce (2018,
2005).

## Acknowledgements

This work was supported by the European Research Council (ERC) under the EU’s
Horizon 2020 research & innovation programme, grant agreement No. 101019375
(Whither Music?), and the Federal State of Upper Austria (LIT AI Lab).

## References

  * Allegraud et al. [2019] Pierre Allegraud, Louis Bigo, Laurent Feisthauer, Mathieu Giraud, Richard Groult, Emmanuel Leguy, and Florence Levé.  Learning sonata form structure on mozart’s string quartets.  Transactions of the International Society for Music Information Retrieval (TISMIR), 2(1):82–96, 2019.
  * Atzmon et al. [2018] Matan Atzmon, Haggai Maron, and Yaron Lipman.  Point convolutional neural networks by extension operators.  ACM Transactions on Graphics, 37(4), 2018.
  * Bresson and Laurent [2017] Xavier Bresson and Thomas Laurent.  Residual gated graph convnets.  In Proceedings of the International Conference on Learning Representations, 2017.
  * Chawla et al. [2002] Nitesh V Chawla, Kevin W Bowyer, Lawrence O Hall, and W Philip Kegelmeyer.  Smote: synthetic minority over-sampling technique.  Journal of artificial intelligence research, 16:321–357, 2002.
  * Chen et al. [2018] Tsung-Ping Chen, Li Su, et al.  Functional Harmony Recognition of Symbolic Music Data with Multi-task Recurrent Neural Networks.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2018.
  * Del Barrio et al. [2018] Eustasio Del Barrio, Juan A Cuesta-Albertos, and Carlos Matrán.  An optimal transportation approach for assessing almost stochastic order.  In The Mathematics of the Uncertain, pages 33–44. Springer, 2018.
  * Deutsch [2013] Diana Deutsch.  Psychology of music.  Elsevier, 2013.
  * Di Giorgi et al. [2021] Bruno Di Giorgi, Matthias Mauch, and Mark Levy.  Downbeat tracking with tempo-invariant convolutional neural networks.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2021.
  * Duane and Pardo [2009] Ben Duane and Bryan Pardo.  Streaming from midi using constraint satisfaction optimization and sequence alignment.  In Proceedings of the International Computer Music Conference (ICMC), 2009.
  * Elowsson and Friberg [2019] Anders Elowsson and Anders Friberg.  Modeling Music Modality with a Key-Class Invariant Pitch Chroma CNN.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2019.
  * Foscarin et al. [2021] Francesco Foscarin, Nicolas Audebert, and Raphaël Fournier-S’Niehotta.  Pkspell: Data-driven pitch spelling and key signature estimation.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2021.
  * Foscarin et al. [2022] Francesco Foscarin, Emmanouil Karystinaios, Silvan David Peter, Carlos Cancino-Chacón, Maarten Grachten, and Gerhard Widmer.  The match file format: Encoding alignments between scores and performances.  In Proceedings of the Music Encoding Conference (MEC), 2022.
  * Fradet et al. [2021] Nathan Fradet, Jean-Pierre Briot, Fabien Chhel, Amal El Fallah Seghrouchni, and Nicolas Gutowski.  MidiTok: A python package for MIDI file tokenization.  In Late-Breaking Demo Session of the International Society for Music Information Retrieval Conference (ISMIR), 2021.
  * Giraud et al. [2015] Mathieu Giraud, Richard Groult, Emmanuel Leguy, and Florence Levé.  Computational Fugue Analysis.  Computer Music Journal, 39(2):77–96, 2015.
  * Hamilton et al. [2017] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  Advances in Neural Information Processing Systems, 30, 2017.
  * Harte [2010] Christopher Harte.  Towards automatic extraction of harmony information from music signals.  PhD thesis, Queen Mary, University of London, 2010.
  * Hentschel et al. [2021] Johannes Hentschel, Markus Neuwirth, and Martin Rohrmeier.  The annotated mozart sonatas: Score, harmony, and cadence.  Transactions of the International Society for Music Information Retrieval (TISMIR), 4(1), May 2021.
  * Hsiao and Su [2021] Yo-Wei Hsiao and Li Su.  Learning note-to-note affinity for voice segregation and melody line identification of symbolic music data.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), pages 285–292, 2021.
  * Hu et al. [2019] Weihua Hu, Bowen Liu, Joseph Gomes, Marinka Zitnik, Percy Liang, Vijay Pande, and Jure Leskovec.  Strategies for pre-training graph neural networks.  In International Conference on Learning Representations (ICLR), 2019.
  * Jeong et al. [2019] Dasaem Jeong, Taegyun Kwon, Yoojin Kim, and Juhan Nam.  Graph neural network for music score data and modeling expressive piano performance.  In Proceedings of the International Conference on Machine Learning (ICML), pages 3060–3070. PMLR, 2019.
  * Karystinaios and Widmer [2022] Emmanouil Karystinaios and Gerhard Widmer.  Cadence detection in symbolic classical music using graph neural networks.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2022.
  * Karystinaios and Widmer [2023] Emmanouil Karystinaios and Gerhard Widmer.  Roman Numeral Analysis with Graph Neural Networks: Onset-wise Predictions from Note-wise Features.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2023.
  * Karystinaios et al. [2023] Emmanouil Karystinaios, Francesco Foscarin, and Gerhard Widmer.  Musical Voice Separation as Link Prediction: Modeling a Musical Perception Task as a Multi-Trajectory Tracking Problem.  In Proceedings of the Joint Conference on Atrificial Intelligence (IJCAI), 2023.
  * Kermarec et al. [2022] Mathieu Kermarec, Louis Bigo, and Mikaela Keller.  Improving tokenization expressiveness with pitch intervals.  In Late-Breaking Demo Session of the International Society for Music Information Retrieval Conference (ISMIR), 2022.
  * Klapuri and Davy [2006] Anssi Klapuri and Manuel Davy.  Signal Processing Methods for Music Transcription.  Springer, 2006.
  * Lattner et al. [2018] Stefan Lattner, Maarten Grachten, and Gerhard Widmer.  Learning Transposition-Invariant Interval Features from Symbolic Music and Audio.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2018.
  * Lattner et al. [2019] Stefan Lattner, Monika Dörfler, and Andreas Arzt.  Learning Complex Basis Functions for Invariant Representations of Audio.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2019.
  * López et al. [2021] Néstor Nápoles López, Mark Gotham, and Ichiro Fujinaga.  AugmentedNet: A Roman Numeral Analysis Network with Synthetic Training Examples and Additional Tonal Tasks.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2021.
  * McLeod and Rohrmeier [2021] Andrew Philip McLeod and Martin Alois Rohrmeier.  A modular system for the harmonic analysis of musical scores using a large vocabulary.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2021.
  * Micchi et al. [2020] Gianluca Micchi, Mark Gotham, and Mathieu Giraud.  Not all roads lead to Rome: Pitch representation and model architecture for automatic harmonic analysis.  Transactions of the International Society for Music Information Retrieval (TISMIR), 3(1):42–54, 2020.
  * Nakamura et al. [2016] Eita Nakamura, Masatoshi Hamanaka, Keiji Hirata, and Kazuyoshi Yoshii.  Tree-structured probabilistic model of monophonic written music based on the generative theory of tonal music.  In Proceedings of the International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 276–280. IEEE, 2016.
  * Pearce [2005] Marcus Pearce.  The Construction and Evaluation of Statistical Models of Melodic Structure in Music Perception and Composition.  PhD thesis, City University of London, UK, 2005.
  * Pearce [2018] Marcus T. Pearce.  Statistical learning and probabilistic prediction in music cognition: Mechanisms of stylistic enculturation.  Annals of the New York Academy of Sciences, 1423(1):378–395, 2018.
  * Satorras et al. [2021] Victor Garcia Satorras, Emiel Hoogeboom, and Max Welling.  E(n) equivariant graph neural networks.  In Proceedings of the International Conference on Machine Learning (ICML), 2021.
  * Schlichtkrull et al. [2018] Michael Sejr Schlichtkrull, Thomas N. Kipf, Peter Bloem, Rianne van den Berg, Ivan Titov, and Max Welling.  Modeling relational data with graph convolutional networks.  In Proceedings of the Semantic Web International Conference, ESWC, volume 10843 of Lecture Notes in Computer Science, pages 593–607. Springer, 2018.
  * Sears et al. [2017] David RW Sears, Andreas Arzt, Harald Frostel, Reinhard Sonnleitner, and Gerhard Widmer.  Modeling harmony with skip-grams.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2017.
  * Wang et al. [2019a] Yue Wang, Yongbin Sun, Ziwei Liu, Sanjay E. Sarma, Michael M. Bronstein, and Justin M. Solomon.  Dynamic graph CNN for learning on point clouds.  ACM Transactions on Graphics, 38(5):146:1–146:12, 2019.
  * Wang et al. [2019b] Yue Wang, Yongbin Sun, Ziwei Liu, Sanjay E Sarma, Michael M Bronstein, and Justin M Solomon.  Dynamic graph cnn for learning on point clouds.  ACM Transactions on Graphics (tog), 38(5):1–12, 2019.
  * Wu et al. [2021] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Philip S. Yu.  A comprehensive survey on graph neural networks.  IEEE Transactions on Neural Networks and Learning Systems, 32(1):4–24, 2021.
  * Zhang et al. [2023] Huan Zhang, Emmanouil Karystinaios, Carlos Cancino-Chacón, Simon Dixon, and Gerhard Widmer.  Symbolic music representations for classification tasks: A systematic evaluation.  In Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2023.
  * Zou et al. [2019] Difan Zou, Ziniu Hu, Yewen Wang, Song Jiang, Yizhou Sun, and Quanquan Gu.  Layer-dependent importance sampling for training deep and large graph convolutional networks.  Advances in neural information processing systems, 32, 2019.