# Investigation of Phonon Lifetimes and Magnon-Phonon Coupling in YIG/GGG
Hybrid Magnonic Systems in the Diffraction Limited Regime

Manoj Settipalli  Ann and H.J. Smead Aerospace Engineering Sciences,
University of Colorado Boulder, Boulder, Colorado 80303, USA  Sanghamitra
Neogi<EMAIL_ADDRESS>Ann and H.J. Smead Aerospace Engineering
Sciences, University of Colorado Boulder, Boulder, CO 80303, USA

###### Abstract

Quantum memories facilitate the storage and retrieval of quantum information
for on-chip and long-distance quantum communications. Thus, they play a
critical role in quantum information processing (QIP) and have diverse
applications ranging from aerospace to medical imaging. It is well established
that quantized vibrations (phonons) in mechanical oscillators can behave
quantum mechanically and play an important role in QIP. Bulk acoustic wave
(BAW) phonons, which vibrate within the bulk of a material, are promising
candidates for storing quantum information due to their long lifetimes. In
this work, we investigate a hybrid magnonic system, where BAW phonons are
excited in a Gadolinium Iron Garnet (GGG) thick film via coupling with
quantized electron spin-waves (magnons) in a Yttrium Iron Garnet (YIG) thin
film. Recent experiments on a millimeter scale YIG/GGG device show that the
memories are limited by the phonon lifetime of 0.2 $\mu$s at room temperature.
The phonon lifetime is expected to be longer in the milliKelvin regime.
However, the phonon lifetime could be limited by diffraction in that regime,
especially for small scale devices. A complete understanding of the
diffraction-limited performance of the hybrid magnonic devices is not
available. We present theoretical and numerical analyses of the diffraction-
limited BAW phonon lifetimes, modeshapes, and their coupling strengths to
magnons in planar and confocal high-overtone bulk aoustic wave resonator
(HBAR) structures. We utilize (1) Fourier beam propagation and (2) Hankel
transform eigenvalue problem approaches and discuss the effectiveness of the
two methods to predict the HBAR BAW phonons. Our analyses predicts the
diffraction-limited phonon lifetime to be on the order of 100 milliseconds in
confocal HBAR structures. We illustrate that a focusing dome could
significantly improve the performance of hybrid magnonic HBAR structures, for
quantum memory and transduction applications.

##  I Introduction

Hybrid quantum systems constitute a rapidly emerging research field due to
their integration of diverse platforms, and ability to synergistically address
individual limitations. Among various systems explored, hybrid magnonic
systems offer unique advantages and as a result, received significant
attention in recent years [1, 2]. In these systems, quantized spin waves
(magnons) coherently couple with other information carriers, such as photons
and phonons. The coupling presents avenues for both fundamental investigations
and practical implementations [3, 4, 5, 6, 7, 8, 9, 10, 11]. Hybrid magnonic
systems often leverage magnetic materials with high spin density, such as
yittrium iron gartnet (YIG). Such materials enable strong coupling of magnons
with microwave photons trapped in a cavity, with coupling strengths surpassing
their respective dissipation rates [3, 4, 5, 6, 7]. Several devices have been
proposed that benefit from the strong magnon-photon coupling, with
applications ranging from microwave-to-optical transduction [12, 13], quantum
magnonics [14, 15, 16], to dark matter detection [17, 18] and more [19, 20,
21, 22, 23, 24].

In addition to magnon-photon coupling, magnon-phonon coupling in hybrid
magnonic devices has attracted attention for coherent and quantum information
processing applications [11, 25, 26, 27, 28]. Magnon-phonon coupling has
already shown success for classical information processing applications [29,
30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41]. Phonons are particularly
attractive due to their long lifetime compared to other information carriers
[42, 43, 44, 45]. Coupling of magnons with mechanical phonons in YIG spheres
has been shown to enable magnon-phonon entanglement [46], nonreciprocal phonon
propagation [47, 48], and other applications. Separately, high-overtone bulk
acoustic wave resonator (HBAR) phonons, which operate in the GHz regime, have
been shown to enable resonant magnon-phonon coupling via the magnetoelastic
effect [49]. The HBAR phonons are commonly supported by structures composed of
YIG thin-films [50, 51, 52, 53] bonded to a thicker substrate, such as a
gadolinium gallium garnet (GGG) film, which can host evenly distributed HBAR
resonances. A recent study demonstrated a triply resonant photon-magnon-phonon
system using a YIG/GGG hybrid magnonic HBAR device [54]. Such resonant
coupling can enable long-lasting multimode phononic quantum memories and
realize quantum information processing and transduction applications. In this
article, we investigate a class of planar and confocal HBAR YIG/GGG hybrid
magnonic structures for quantum transduction applications.

The performance figure of merit for magnon-phonon transduction can be measured
using cooperativity, $C=4g_{mb}^{2}/\kappa_{m}\kappa_{b}$, where $g_{mb}$,
$\kappa_{m}$, and $\kappa_{b}$ are magnon-phonon coupling strength, magnon
dissipation rate, and phonon dissipation rate, respectively. The lifetime of
magnons and phonons are given by $\tau_{m}=1/\kappa_{m}$ and
$\tau_{b}=1/\kappa_{b}$, respectively. It is commonly accepted that the phonon
lifetimes surpass the magnon lifetimes in hybrid magnonic devices. Past
studies reported the magnon and phonon lifetime in YIG/GGG HBAR systems to be
$\tau_{m}=0.07\mu$s and $\tau_{b}~{}\sim 0.25\mu$s, respectively [53, 54]. It
is desirable to improve the system lifetime much beyond $\sim 0.25\mu$s for
quantum memories and other information processing applications. At given
temperature and frequency operating conditions, the phonon lifetime could be
limited by material and diffraction losses. However, a complete understanding
about the loss mechanisms in devices with given geometry or diverse operating
conditions is not available. Phonon lifetime at room temperature is primarily
limited by acoustic attenuation due to phonon-phonon interactions [55, 56].
The acoustic attenuation effects, $\alpha\propto 1/\tau_{b}$, has been shown
to follow a $T^{4}$ temperature dependence. Due to the $T^{4}$ behavior, these
effects are less significant at low temperatures, such as mK regime where many
qubit systems operate. Ideally, $\tau_{b}$ can increase by multiple orders of
magnitude at cryogenic temperatures. However, it is important to consider the
diffraction losses to decide if low-T operations are a viable solution to
generate long-lifetime phonons. In this study, we investigate the effects of
diffraction losses on HBAR phonon modes in hybrid magnonic structures and
leave the study of material losses for future research. One of our main
objectives is to establish a predictive approach to analyze diffraction losses
in HBAR structures. However, we like to emphasis that the knowledge of
material losses in HBAR structures across various temperature and frequency
operating conditions is limited. There is a strong need for experimental
characterization of acoustic attenuation in YIG/GGG material systems to unlock
their full potential.

One important design aspect of hybrid magnonic devices is its scalability: the
number of memory or transduction centers that could be incorporated into a
single chip, in addition to other on-chip circuitry. The YIG film represents
the transduction center in the YIG/GGG HBAR structures of interest. The
smaller the lateral area of the YIG film, the greater the number of such
centers in a single chip and the scalability. We can define the scalability
factor, $S$, as $S=\frac{A^{0}_{d}}{A_{d}}$. Here, $A_{d}$ is the lateral area
of the YIG film for planar HBAR structures, and the cross-sectional area of
the confocal dome surface for CHBAR structures in our study, respectively.
$A^{0}_{d}$ represents a reference value for the device. We consider the
smallest known lateral area reported in literature to be the reference,
$A^{0}_{d}=0.8\times 0.9$ mm${}^{2}=0.72$ mm2 [54]. The scalability of the
device can be increased by reducing the lateral area of the YIG film. However,
the reduced aperture will lead to high diffraction and affect phonon lifetime.
In addition, a reduced magnon-phonon overlap in the YIG region could lead to a
weaker magnon-phonon coupling. The reduced lifetime and the coupling will
decrease the cooperativity. Therefore, a complete understanding of the trade-
off between scalability and diffraction-limited performance is essential to
obtain long-living scalable phononic quantum memories using HBAR structures.
Note that we only consider the lateral scalability and use the same thickness,
$527.2\mu$m, for all structures, to keep the phonon free spectral range
relatively unchanged.

In this study, we investigate HBAR phonon modes, their diffraction-limited
lifetimes, and the magnon-phonon coupling strengths in planar and confocal
HBAR YIG/GGG hybrid magnonic structures and identify strategies to improve the
performance and scalability of these structures. We use the magnon lifetime to
be $0.07\mu$s for all structures considered in this study [54]. The phonon
lifetime $\tau_{b}$ or the quality factor ($Q=\omega\tau_{b}$) is of
particular interest due to its direct relevance to the quantum memory storage
time. Henceforth, we drop the subscript $b$ from $\tau_{b}$, unless otherwise
needed to distinguish it from magnon and photon lifetimes. We consider the
magnon mode to be the fundamental mode or the Kittel mode, which has a well-
degined analytical function for circular discs. However, no such analytical
description exists for the HBAR phonon modes. We evaluate the existing
numerical methods and discuss the appropriate methods to reliably model the
HBAR phonons modes. One method estimated phonon lifetimes in planar HBARs by
decomposing the initial beam in a Bessel function basis. They calculated the
overlap of the propagating Bessel functions with the initial profile after
every round trip and extracted the lifetime from an exponential fit [57].
However, the predicted lifetime was smaller than the experimental observations
because they ignored the effects of the lateral confinement due to the
transducer.

Finite element analysis (FEA), as implemented in the COMSOL Multiphysics [58]
software, has been a popular method to study HBAR phonons in planar and
confocal geometries. A recent study estimated the diffraction loss in an
epitaxial planar HBAR structure by measuring the power received at the
opposite end to that of the actuator [59]. However, they considered half a
round trip and also ignored the effect of localization due to the transducer.
In general, acoustic wave amplitudes decay non-exponentially with time due to
the mode mismatch between the input acoustic beam profile and the phonon mode
of interest. Hence, it is important to consider how transducer localization
affects over multiple round trips to accurately determine the characteristic
decay of the mode of interest, instead of only half a round trip. Another
study used FEA to predict the modeshapes of a CHBAR for qubit coupling
applications [60]. They considered phonon modes at low-overtones or long
wavelengths, hence it was possible to perform 3D FEA simulations for this
system. However, as we will discuss later, it becomes intractable to simulate
bulk acoustic waves with overtone numbers $n{\sim}3000$ using FEM since
sampling required for these small wavelength modes is on the order of
$~{}{\sim}$ billion nodes. All these studies discussed longitudinal HBAR
phonons. However, forward volume magnetostatic modes, such as Kittel modes,
strongly couple to shear phonons [53, 54, 61]. The FEA analysis becomes
particularly expensive for shear phonons since they are not axi-symmetric. One
cannot leverage the axi-symmetric 2D FEA modelling available in COMSOL
Multiphysics for these systems. A past study simulated a 2D YIG/GGG HBAR
system in COMSOL without leveraging axi-symmetry [54]. They ignored the
different impacts of Gouy phase effects between 2D and 3D diffraction [62]. It
is necessary to consider such effects to accurately predict lifetimes and mode
shapes.

Due to these reasons, we explore other approaches to model shear phonon modes.
We consider the Fourier beam propagation method (FBPM) [63, 57] which follows
a Fox-Li-like [64] iterative approach to obtain the phonon modes and works
with a plane-wave basis set. Although past studies used FBPM for longitudinal
phonon modes, we illustrate that it can effectively model the shear phonon
modes of interest. We provide s detailed description of this method and
discuss an adaptive algorithm that allows to overcome some of the challenges
of the standard FBPM method. In addition we consider another method that uses
Hankel (HK) transform, which is a Fourier transfrom for axi-symmetric systems,
and works with a Bessel function basis. Thus far, HK method has only been
implemented for Fabry-Perot optical cavities [65, 66]. We adapt this approach
to simulate YIG/GGG HBAR structures by leveraging their axi-symmetry and
isotropic material properties.

##  II Methods

###  II.1 HBAR Configurations

Figure 1 shows the two representative HBAR structures investigated in this
work. The structure in Figure 1(a) includes a thick GGG film joined together
with a YIG thin film at the bottom. We refer to this structure as the planar
HBAR structure since it has planar top and bottom surfaces. Figure 1(b) is a
focusing HBAR structure that includes a GGG dome structure at the top and a
planar bottom surface. We refer to this structure as the confocal HBAR (CHBAR)
structure. All configurations considered here have a GGG layer of thickness,
$t_{\text{GGG}}=527\mu$m. This value refers to the length from the bottom
surface of the GGG film to the mid point of the dome for the confocal HBAR, as
shown. All YIG films are circular films or discs with thickness,
$t_{\text{YIG}}=200$ nm, and radius of cross-section, $R_{\text{YIG}}$. Thus,
the total HBAR device thickness is
$t_{\text{HBAR}}=t_{\text{GGG}}+t_{\text{YIG}}=527.2\mu$m. We consider planar
and confocal HBAR structures with several $R_{\text{YIG}}$’s ranging from
10$\mu$m to 200$\mu$m. The width $W$ of the structure is 1200$\mu$m unless
otherwise mentioned. The radius of cross-section of the dome,
$R_{\text{cross}}$, and its radius of curvature, $R_{\text{curv}}$, vary for
different structures considered. Altogether, we consider 8 planar HBARs with
varying $R_{\text{YIG}}$, 18 CHBARs with varying $R_{\text{curv}}$ and fixed
$R_{\text{cross}}$, and 10 CHBARs with varying $R_{\text{cross}}$ and fixed
$R_{\text{curv}}$.

Figure 1: Representative HBAR configurations: (a) Planar HBAR structure
composed of a GGG thick film and a YIG thin film. (b) Confocal HBAR structure
with a dome surface at the top and a planar bottom surface. $x$ axis is along
the out-of-plane direction while $y$ and $z$ axes point along the lateral and
the normal direction of the structures, respectively. The origin of the
coordinate axes is at the center of the bottom YIG surface.

###  II.2 Magnon Modes

The YIG thin film hosts magnons, generated by a static, $\textbf{B}_{0}$, and
an oscillating RF magnetic field, $\textbf{B}_{\textbf{RF}}$. $\textbf{B}_{0}$
is applied along the normal direction of the YIG film, while the oscillating
$\textbf{B}_{\textbf{RF}}$ acts in the YIG plane. The magnetic fields,
$\textbf{B}_{\textbf{0}}$ and $\textbf{B}_{\textbf{RF}}$, generate forward
volume magnetostatic modes, as shown in Fig. 1. These modes precess in the
$x-y$ plane. The dynamic magnetization in the YIG region is given by
$\textbf{m}_{\text{YIG}}(x,y)=m_{0}(x,y)(\text{cos}(\omega
t)\textbf{i}+\text{sin}(\omega t)\textbf{j})$. Here, $m_{0}(x,y)$ represents
the magnon modeshape and $\omega$ is the frequency of precession. We practice
the following convention throughout the article, bold-font symbols represent
vector fields and corresponding normal-font symbols represent scalar values,
respectively. We expect un-pinned spin waves at the top and bottom surfaces.
This prediction can be justified using the following two arguments: (1)
exchange interactions are expected to dominate at a thickness on the order of
200 nm and below and (2) as a consequence, the pinning effect is expected to
disappear below a critical width on the order of 200 nm [67]. However, the
magnetic dipolar interactions are expected to dominate at length scales on the
order of micrometers, over the exchange interactions [67]. As a result, we
expect full-pinning of the magnetic spin waves at the lateral boundaries of
the YIG discs of radii ranging from 10$\mu$m to 200$\mu$m. Considering these
two expected results, we assume that $\textbf{m}_{\text{YIG}}$ is constant
throughout the thickness of the YIG film along the $z$ axis. Accordingly, we
describe the modeshape of the pinned spin waves, $m_{0}(x,y)$, using truncated
Bessel functions [68, 69]:

$m_{0}(x,y)=\begin{cases}J_{0}\left(\frac{R}{R_{\text{YIG}}}\zeta_{j}\right),&\text{if
}R\leq R_{\text{YIG}}\\\ 0,&\text{otherwise}\end{cases}$ (1)

where radial coordinate $R=\sqrt{x^{2}+y^{2}}$ and $J_{0}$ is the
$0^{\text{th}}$ order bessel function. The zeroes of $\zeta_{j}$ with
$j={0,1,2,...}$, correspond to the fundamental and higher-order magnon modes
respectively. In this article, we focus on the fundamental magnon mode or the
Kittel mode, whose amplitude is represented by the truncated function,
$J_{0}(\frac{R}{R_{\text{YIG}}}\zeta_{0})$. The modeshape is constant along
the $z$-direction.

###  II.3 Phonon Modes

The precessing forward volume magnon modes, $\textbf{m}_{\text{YIG}}(x,y)$,
generate precessing shear deformations in the YIG region. The dynamic shear
strain results in a circularly polarized chiral phonon traveling wave in the
GGG region. The helicity of the chiral phonon mode is determined by the
precession direction of the Kittel mode. The dynamic shear displacements can
be expressed as $\textbf{u}_{\text{YIG}}(x,y)=u_{0}(x,y)(\text{cos}(\omega
t)\textbf{i}+\text{sin}(\omega t)\textbf{j})$, where $u_{0}(x,y)$ is the
phonon modeshape and $\omega$ is the precession frequency. In the following,
we refer to the shear displacements, $\textbf{u}_{\text{YIG}}(x,y)$, as
$\textbf{u}_{0}(x,y)$. The traveling wave undergoes reflections at the top GGG
surface and the reflected waves interfere with the forward traveling wave.
When the forward and reflected helical propagating waves interfere, they form
rotating standing shear wave modes at certain frequencies. While the traveling
chiral phonons are helical, the standing waves are not helical. The top GGG
surface does not induce a $\pi$ phase shift to the reflected wave, unlike
circularly polarized light reflecting off a mirror. As a result, we obtain
standing shear wave modes with zero net helicity. In this article, we use (1)
Fourier beam propagation and (2) Hankel transform eigenvalue problem
approaches to analyze the shear wave phonon modes in the chosen HBAR
configurations.

####  II.3.1 (1) Phonon Modeshape Analysis: Fourier Beam Propagation

The Fourier beam propagation method (FBPM), also known as the angular spectrum
method, predicts the field displacements or profiles of propagating waves at
distances away from the aperture [63, 70]. It has been widely used to analyze
beam propagation in the field of optics [70]. More recently, the FBPM is used
to study HBAR phonons in planar [57] and CHBAR structures [63], adapting an
iterative method similar to the Fox-Li approach [64]. The advantage of FBPM
lies in its simplicity and the ability to predict the field profiles at any
target distances from the source without needing to calculate the behavior at
intermediate distances. The mathematical formulation of FBPM for acoustic
waves in anisotropic systems is well established [63]. Here, we implement a
reformulated approach in which the propagation is calculated using projector
and propagator operators. This reformulation allows us to achieve a seven-fold
speed up of the computation time. Our implementation is particularly
advantageous for isotropic systems such as YIG and GGG.

The iterative procedure begins with an initial input displacement field which
we assume to be

$u_{0}^{(1)}(x,y)=\begin{cases}J_{0}\left(\frac{R}{R_{\text{YIG}}}\zeta_{0}\right),&\text{if
}R\leq R_{\text{YIG}}\\\ 0,&\text{otherwise}\end{cases}$ (2)

where the superscript represents iteration index ($i=1$) and the subscript in
$u^{(1)}_{0}(x,y)$ refers to $z=0$. Note that the lateral phonon modeshape is
the same as the Kittel modeshape mentioned earlier in Eq. 1. In FBPM, the
input beams for initial and the subsequent iterations are decomposed into
plane-waves. The field displacements at intermediate distances is then
obtained by multiplying phase factors to the decomposed beam. We decompose the
input acoustic beam into plane-waves using Fourier decomposition:

$\tilde{\textbf{u}}^{(i)}_{0}(k_{x},k_{y})=\text{FFT}[\textbf{u}^{(i)}_{0}(x,y)].$
(3)

Here, $\tilde{\textbf{u}}^{(i)}_{0}$ is defined on a $N\times N$
$(k_{x},k_{y})$ grid which is the Fourier conjugate of the spatial $N\times N$
$(x,y)$ grid. The index $i$ represents the iteration number. We calculate the
projections of $\tilde{\textbf{u}}^{(i)}_{0}(k_{x},k_{y})$ along the different
polarization directions to obtain the amplitudes, $A^{(i)}_{m}$, of the shear
($m=1,2$) and the longitudinal ($m=3$) modes:

$\displaystyle A^{(i)}_{m}=\textbf{P}_{m}\cdot\tilde{\textbf{u}}^{(i)}_{0},$
(4a) with
$\displaystyle\textbf{P}_{m}=\frac{\Sigma_{cyc}(-1)^{|\text{sgn}(i-m)|}\hat{d}_{i}(\hat{d}_{m}\cdot\hat{d}_{i}-\hat{d}_{j}\cdot\hat{d}_{k})}{1+2\Pi_{cyc}(\hat{d}_{i}\cdot\hat{d}_{j})-\Sigma_{cyc}(\hat{d}_{i}\cdot\hat{d}_{j})^{2}},$
(4b)

where $\textbf{P}_{m}$ is the projector onto polarization direction $m$. Here,
$\hat{d}_{m}$ is the polarization vector for the plane-waves propagating along
$(k_{x},k_{y},k_{z,m})$, $(i,j,k)$ refers to the Cartesian directions,
$(i,j,k)=(1,2,3)$ and $\Sigma_{cyc}$ and $\Pi_{cyc}$ are cyclic sum and
product operators, respectively, that cycle through variables $i$, $j$, and
$k$. We use the amplitudes $A^{(i)}_{m}(k_{x},k_{y})$, to obtain the
displacement field of the propagated beam, starting from the input beam,
$u^{(i)}_{0}(x,y)$. For a wave originating in the YIG thin film, the
propagation distance to reach the upper GGG surface of the HBAR structures is
$t_{\text{HBAR}}$, as shown in Fig. 1. The displacement field of the
propagated beam at the upper GGG surface is given by

$\displaystyle\textbf{u}^{(i)}_{t_{\text{HBAR}}}(x,y)=\text{IFFT}[\Sigma_{m}A^{(i)}_{m}G_{m}],$
(5a) with $\displaystyle
G_{m}(k_{x},k_{y})=\hat{d}_{m}(k_{x},k_{y})e^{ik_{z,m}(k_{x},k_{y})t_{\text{HBAR}}}.$
(5b)

Here, $G_{m}(k_{x},k_{y})$ is the propagator for plane-waves traveling along
$(k_{x},k_{y},k_{z,m})$ with polarization $m$. Typically, $k_{z,m}$ is derived
from their respective slowness surfaces [63] for each polarization and values
of $(k_{x},k_{y})$. However, the functional relation can be simplified to
$k_{z,m}=\sqrt{\frac{\omega^{2}}{v^{2}_{m}}-k^{2}_{x}-k^{2}_{y}}$, for
isotropic dispersions. We use this relationship for the propagating waves in
our isotropic YIG/GGG HBAR structures. Here, $\omega$ is the frequency of the
initial wave and $v_{m}$ is the velocity of phonons with polarization $m$.
Note that we need to compute $k_{z,m}$, $A_{m}$ (Eq. 4), and
$\textbf{u}_{t_{\text{HBAR}}}(x,y)$ (Eq. 5) only during the first round trip,
for a given $\omega$. This aspect results in the seven-fold computational
speed-up mentioned earlier. We show the properties of YIG and GGG in Table. 1.
It can be noted that the mechanical properties of YIG and GGG are similar. Due
to this reason, we use GGG values for both YIG and GGG, for simplicity. This
approximation can be further justified by considering that our structures are
composed of GGG thick films with ultrathin YIG films bonded to it.

Table 1: Material properties of YIG and GGG. | Young’s modulus | Poisson’s ratio | Density
---|---|---|---
| E (Pa) | $\nu$ |  $\rho$ (kg/m3)
YIG [71] | 0.2 $\times$ 1012 | 0.29 | 5170
GGG [72] | 0.222 $\times$ 1012 | 0.28 | 7080

We choose the longitudinal polarization $\hat{d}_{3}$ to be along the unit
vector k and the shear polarizations, $\hat{d}_{1,2}$, to be two mutually
perpendicular unit vectors perpendicular to k. Note that for isotropic systems
with mutually perpendicular $d_{m}$s, like what is considered in this article,
$A_{m}$ reduces to $A_{m}=\hat{d}_{m}\cdot\tilde{\textbf{u}}$ at each
$(k_{x},k_{y})$.

The propagated plane-waves are periodic in the direction transverse to the
beam propagation direction. When the diffracted waves reach the transverse
boundaries of the computational domain, they introduce undesired reflections.
To avoid these reflections, one needs to consider a sufficiently wide
computational domain that can contain the waves even after multiple
reflections. However, such large domains can significantly increase
computational costs. An alternative approach is to introduce absorbing
boundary regions that completely attenuate any waves that enter these regions
[57]. In this study, we implement absorbing boundaries of thickness
$W_{\text{Abs}}=50\times\lambda$ (wavelength of the input field). We show the
absorbing boundaries in in Fig. 1, using blue shaded regions. To simulate the
effect of absorbing boundaries, we multiply
$\textbf{u}^{(i)}_{t_{\text{HBAR}}}$ (Eq. 5) with a reflection operator
$R_{t_{\text{HBAR}}}$, defined as

$R_{t_{\text{HBAR}}}(x,y)=\begin{cases}1,&\text{if }R\leq W_{\text{eff}}/2\\\
0,&\text{otherwise}\end{cases}$ (6)

where $W_{\text{eff}}=W-2W_{\text{Abs}}$ is the effective width of the
simulation window without the absorbing boundaries. We propagate the
attenuated reflected wave further through a distance $t_{\text{HBAR}}$ to
complete a full round trip. We implement the beam propagation by following the
approach outlined in Eqs. 3 \- 5b. After every round trip, we multiply the
resulting complex displacement field by an additional phase $R_{0,m}$ for each
polarization $m$:

$R_{0,m}(x,y)=\begin{cases}e^{i2k_{z0,m}t_{\text{YIG}}},&\text{if
}R_{\text{YIG}}\leq R\leq W_{\text{eff}}/2\\\ 1,&\text{if }R\leq
R_{\text{YIG}}\\\ 0.&\text{otherwise}\end{cases}$ (7)

where $k_{z0,m}=k_{z,m}(k_{x}=0,k_{y}=0)$. The phase factor is introduced due
to the finite width of the YIG film compared to the GGG width, $W$. However,
it is worth mentioning that this approximation is more appropriate for low-
diffraction cases. We used this for all cases considered to simplify the
analysis. We use the resulting displacement field as the new input field for
the next round trip iteration. We repeat this process for $N$ round trips. We
calculate the complex sum of the displacement fields,
$\textbf{U}_{0}(x,y)=\Sigma^{N}_{i=1}\textbf{u}^{(i)}_{0}(x,y)$ at $z=0$, to
include the effect of interferences. We use the interference sum
$\textbf{U}_{0}(x,y)$ obtained after $N$ round trips as the initial beam of
the next restart, $\textbf{u}^{(1)}_{0}(x,y)=\textbf{U}_{0}(x,y)$, and restart
the process. Such a restart process using interference sum as an input ensures
fast convergence to the desired mode. The other mode components are attenuated
in the interference sum through destructive interference. We continue this
process until the varation of the modeshape is within a chosen tolerance. We
obtain a converged standing wave with displacement field,
$\textbf{u}_{0}(x,y)$, as a final outcome.

The frequency overtones for plane-waves traveling along z direction in the
HBAR structures are expected to be
$\omega_{m,n}=2\pi\times\frac{nv_{m}}{2t_{\text{HBAR}}}$, with
$n=1,2,3,...\infty$. The thickness of the structure is $t_{\text{HBAR}}$,
velocity of wave is $v_{m}$ and $m$ represents polarization. The overtones are
separated by $\frac{v_{m}}{2t_{\text{HBAR}}}$, known as the free spectral
range (FSR). However, the diffracting waves traveling in the $z$ direction do
not have well-defined $k_{z,m}$’s leading to Gouy phase effects [62], unlike
plane-waves. The diffraction results in the shift of resonance frequencies
from the monochromatic plane-wave overtones. In this work, we are interested
in identifying the frequencies of the shear modes with $m=2$. The reason is
that the shear modes could couple with the Kittel magnon modes generated in
the YIG thin-film (Eq. 1), in a rotating coordinate system. Note that it
suffices to consider only $m=1$ or $m=2$ as the dominant shear polarization in
the rotating coordinate system. To identify the resonance frequencies, we
sweep through a frequency window of $\omega_{0}\pm 5\times\text{FSR}$ divided
into 50 steps. Here, the frequency of interest, $\omega_{0}=2\pi\times 9.825$
GHz, corresponds to the frequency of the 2960th overtone of a standing plane-
wave in the HBAR structure. We choose this overtone to maintain the same
frequency operating conditions as was considered in a previous article that
investigated similar HBAR structures [54]. For each frequency, we propagate
the beam in the structure and calculate the intensities of the interference
sums ($I=\int_{A}\text{Re}[\textbf{U}_{0}(x,y)]^{2}dA$) at $z=0$ after $N$
round trips. We obtain the resonance frequencies by noting the frequencies for
which the intensities are maximum. We consider a YIG/GGG HBAR structure with
$t_{\text{YIG}}=200$ nm, $R_{\text{YIG}}=200\mu$m, $t_{\text{GGG}}=527\mu$m,
and $L_{x,\text{GGG}}=L_{y,\text{GGG}}=1200\mu$m, for this analysis.

Figure 2: Resonant modes in planar HBAR structures: (a) Multimode phonons,
represented by the overtones of the fundamental mode, and their free spectral
range (FSR). (b) Isometric view of the y-component of $\textbf{U}_{0}(x,y)$,
$\text{U}_{0y}(x,y)$. (c). Top view of $\text{U}_{0y}(x,y)$ showing
localization effect caused by the YIG film. (d) Weighted deviation (WD)
between mode profiles before first and second restarts.

We show the intensities of the propagated beam in Fig. 2(a), for the frequency
range considered. The central intensity peak corresponds to
$\omega_{0}=2\pi\times\frac{2960v_{2}}{2t_{\text{HBAR}}}$, the exact value of
the $2960^{\text{th}}$ shear plane-wave overtone. The peaks form at resonance
frequencies separated by FSR $=3.32$ MHz, indicating the presence of multiple
shear wave overtones around the frequency, $\omega_{0}/2\pi\sim 9.825$ GHz. In
Figs. 2 (b) and (c), we show the $y$-component of the normalized resonant
modeshape, Re[U${}_{0y}(x,y)$] near $\omega_{0}$. We obtain these modeshapes
after 800 round trips. This simulation included one restart after 400 round
trips. This number is arbitrary but we chose after performing multiple tests.
The basis of selection is that the number is large enough such that the
effects from Gouy phase are not significant on the interference sums. We
compute the weighted deviation (WD) between the normalized displacement
fields, obtained before the 1st and the 2nd restarts,
$\textbf{U}^{1}_{0}(x,y)$ and $\textbf{U}^{2}_{0}(x,y)$, respectively:
$|(\textbf{U}^{1}_{0}(x,y)-\textbf{U}^{2}_{0}(x,y))|/|\Sigma_{x,y}\textbf{U}^{1}_{0}(x,y)|$.
We monitor the WD between restarts to check for convergence. Figure 2 (d)
shows that the weighted deviation between the displacement fields. As can be
noted, WD is much less than $10^{-7}$, indicating that the Figs. 2 (b) and (c)
represent displacement fields converged within sufficient numerical tolerance.
Some fringes remain in the modeshape profiles that are possible artifacts of
our numerical analysis. These are the high frequency components resulting from
the sampling of the sharp phase change induced by the YIG film. They are
significantly lower in values, however, could be further reduced by either
low-pass filtering or using a smoother phase change factor than $R_{0,m}(x,y)$
(Eq. 7) used in this work.

However, to obtain converged modes and identify the true modal frequency near
an intensity peak of interest, we further narrow the frequency sweeping range
with finer sampling ($\sim 1$ kHz). To identify the true modal frequency near
$\omega_{0}$, we set the frequency to be
$\omega^{\text{HBAR}}_{0}=\omega_{0}+\delta\omega$, with
$\delta\omega=2\pi\times 2.157$ kHz and restart the iterative process with Eq.
2 as the input beam. If such narrowing is not done, the modeshape can change
significantly after each restart. This is due to the Gouy phase effects the
beam incurs due to diffraction [62], as discussed earlier. This results in the
detuning of overtones with respect to the plane-wave overtones. We illustrate
the detuning effect in Fig. 3. We show the real and imaginary parts of the
complex sums U${}^{(1)}_{0y}(x,y)$ and U${}^{(2)}_{0y}(x,y)$ computed before
first and second restart (after 400 and 800 round trips), respectively. We
calculate the complex sums at $\omega_{0}$ without performing the
$\delta\omega$ correction discussed above. It can be noted that both the real
and imaginary parts of $\text{U}^{(1)}_{0y}$ and $\text{U}^{(2)}_{0y}$ vary
significantly during the round trips. Therefore, this approach does not result
in converged mode shapes.

Figure 3: Lateral mode profiles obtained with FBPM showing Gouy phase effects:
The real (solid) and imaginary (dashed) components of the y-component of the
interference sum before the first $\text{U}^{(1)}_{0y}$ and second restarts
$\text{U}^{(2)}_{0y}$ at
$\omega_{0}=2\pi\times\frac{2960v_{2}}{2t_{\text{HBAR}}}$. The components
deviate significantly between restarts.

.

Although the narrowed frequency sweeping helps with identifying the true
resonant frequency, it is a tedious process to identify the necessary fine
sampling. We need to restart the iterative process multiple times,
particularly for high-diffraction cases. The results shown in Fig. 2,
represent propagating waves in a low diffraction regime with Fresnel number
$N_{F}=213$. However, this work also considers wave propagation in HBAR
structures with $R_{\text{YIG}}$ as low as $10\mu$m corresponding to a very
low Fresnel number, $N_{F}=1.06$ indicating a high-diffraction regime. We
investigate these structures to discuss the effects of reduced actuator
lateral area on the phonon modes and their lifetimes. The reduced area
promises to increase scalability of the devices for memory applications. The
high diffraction regime of these structures introduces significant Gouy phase
effects on the propagating waves and detuning of frequencies. It becomes
challenging to identify a resonant frequency by merely narrowing the frequency
sweeping range. We implement an adaptive version of the FBPM algorithm to
circumvent this challenge. We like to point out that the phonon modeshapes can
be well predicted by using absorbing boundaries in the FBPM technique.
However, their lifetime can be dependent on how the shape and size of the
absorbing boundaries are defined. We do not carry out extensive investigation
to narrow down the choice of the boundaries and left for a future
investigation.

####  II.3.2 Phonon Modeshape: Adaptive Fourier Beam Propagation

The FBPM approach described above can be formulated as an eigenvalue problem
as given below:

$\textbf{RT}\textbf{u}_{\text{0}}(x,y)=\Lambda\textbf{u}_{\text{0}}(x,y)$ (8)

where RT is the round trip operator that includes all the transformations
described in the previous section (Eqs. 3\- 7) and $\Lambda$ is the eigenvalue
corresponding to the eigenvector $\textbf{u}_{0}$. $\Lambda$ is real for a
standing wave mode. A complex $\Lambda$ applies a global phase to the input
beam profile after a round trip, which implies a traveling wave mode. The
phase prevents the waves to interfere constructively over multiple round
trips. It is possible to solve the eigenvalue problem for two-dimensional (2D)
problems (e.g., if HBARs are represented by planar 2D structure), however, it
cannot be directly solved for three-dimensional problem of our interest. We
develop an iterative method to obtain the mode that satisfies Eq. 8 with a
real $\Lambda$. In the adaptive FBPM (a-FBPM) algorithm, we iteratively adjust
the frequencies based on the phase difference incurred over a round trip to
arrive at the desired standing wave modes. We outline the steps of the
iterative method below.

  1. 1.

Start with an input beam, $\textbf{u}^{(i)}_{\text{0}}(x,y)$, that travels in
the HBAR structure, for the $i^{th}$ round trip. The initial input beam
$(i=1)$ has an initial wavelength, $\lambda^{(1)}$ $=\frac{2t}{n}$ and
frequency $\omega^{(1)}=2\pi\times\frac{nv_{2}}{2t}$. Here, $n$ is the
overtone number, $t=t_{\text{HBAR}}$ and $v_{2}$ is the velocity of shear
phonons. The wavelength and frequency are estimated based on how overtones are
estimated for an open-ended column for $n^{\textbf{{th}}}$ overtone. The
initial input modeshape is chosen to be a truncated Bessel function as shown
in Eq. 2.

  2. 2.

Calculate
$\textbf{u}^{(i+1)}_{\text{0}}(x,y)=\textbf{RT}\textbf{u}^{(i)}_{\text{0}}(x,y)$.

  3. 3.

Estimate the real eigenvalue,
$\Lambda^{(i)}=\sqrt{\frac{I^{(i+1)}}{I^{(i)}}}$, where
$I^{(i)}=\int|\textbf{u}^{(i)}_{\text{0}}|^{2}dxdy$ and
$I^{(i+1)}=\int|\textbf{u}^{(i+1)}_{\text{0}}|^{2}dxdy$.

  4. 4.

Calculate the residual,
$\text{R}^{(i)}=||(\textbf{u}^{(i+1)}_{\text{0}}(x,y)-\Lambda^{(i)}\textbf{u}^{(i)}_{\text{0}}(x,y))||/||\Lambda^{(i)}\textbf{u}^{(i)}_{\text{0}}(x,y)||$.

  5. 5.

If R(i) $\leq\text{tol}$, end simulation and output final eigenvalue and
modeshape of resonant modes, else continue to the next step. We choose the
tolerance to be $\text{tol}=10^{-6}$ for all the a-FBPM calculations of this
study.

  6. 6.

If R${}^{(i)}>\text{tol}$, estimate the global phase factor induced by the
round trip operation RT,
$\theta^{(i)}=\text{Arg}(\textbf{u}^{(i)}_{\text{0}}(0,0))-\text{Arg}(\textbf{u}^{(i+1)}_{\text{0}}(0,0))$.

  7. 7.

Set $\omega^{(i+1)}=\omega^{(i)}+\frac{v_{2}\theta^{(i)}}{2t_{\text{HBAR}}}$
and $\lambda^{(i+1)}$ accordingly ($\lambda^{(i+1)}=2\pi
nv_{2}/\omega^{(i+1)}$). The steps 6 and 7 makes the algorithm adaptive.

  8. 8.

Repeat the process until step 5 is satisfied or the maximum number of round
trips $N$ is reached, at which point we set
$\textbf{u}^{(1)}_{0}(x,y)=\textbf{U}_{0}(x,y)=\Sigma^{N}_{n=1}\textbf{u}^{(i)}_{0}(x,y)$,
and we restart the iterative procedure.

Figure 4: Obtaining resonant modes using the adaptive FBPM algorithm: Decay of
residual (R) of modeshapes as a function of number of round trips considered.
Results are shown for two planar HBAR structures with $R_{\text{YIG}}=$ 200
$\mu$m (red, left and bottom axes) and 10 $\mu$m (blue, right and top axes).
The circles with arrows point to the axis corresponding to each plot.

We demonstrate the effectiveness of the adaptive FBPM algorithm in predicting
the resonant modes in HBAR structures in Fig. 4. We show the decay of the
residual, R, as a function of the number of steps. We show the results for two
planar HBAR structures with $R_{\text{YIG}}=$ 200 $\mu$m (red) and 10 $\mu$m
(blue), respectively. We calculate the complex interference sum after every
$NR$ round trips and restart the iterative procedure. The choice of the
number, $NR$, is an informed guess. We choose $NR=400$ and 40 for the HBARs
with $R_{\text{YIG}}=$ 200 $\mu$m and 10 $\mu$m, respectively. We choose a
smaller number of round trips, $NR$, for the 10 $\mu$m case, since the HBAR
with $R_{\text{YIG}}=$ 10 $\mu$m represents a high-diffraction structure
compared to the 200 $\mu$m case and the modeshape decays rapidly for this
case. As shown in Fig. 4, it takes a total of $\sim$650 and $\sim$1200 round
trips (including restarts) to obtain converged resonant modes for the 10
$\mu$m and the 200 $\mu$m case, respectively. The step-like features of the
residual decay coincide with the steps when we compute the interference sums
and use it as input for the next restart of the iterative procedure. The jumps
occur because some errors get canceled due to destructive interference when an
interference sum is computed. In some cases, R can have a momentary rise due
to accumulation of errors from previous round trips, however, the overall
trend continues to decrease ultimately reaching convergence. As shown in Fig.
4, we continue the iterative process till R $<10^{-8}$ to show continued
convergence even after the tolerance is reached.

####  II.3.3 (2) Phonon Modeshape Analysis: Hankel Transform Eigenvalue
Problem

As we discussed previously, the challenges with the standard FBPM approach are
that it is tedious to identify the converged mode by sweeping through
frequencies and obtaining the convergent mode is a slow process. The a-FBPM
approach avoids frequency sweeping and mostly overcomes the slow convergence
issues. However, as we discuss later, a-FBPM still suffers from convergence
challenges for some cases, particularly for the confocal HBAR. Here, we
discuss an alternate method to the FBPM/a-FBPM methods, which leverages the
axi-symmetry of the YIG/GGG HBAR system due to their isotropic material
properties and the circular cross-section of the YIG film. Note that although
the vector field of the shear acoustic phonon mode does not obey axi-symmetry,
the scalar field describing the dominant shear-component of the
$\textbf{u}_{0}$ is expected to obey axi-symmetry. Such a system can be
modelled using the axi-symmetric equivalent of a Fourier transform, the Hankel
(HK) transform, reducing the problem from a from the three-dimensional (3D) to
a two-dimensional (2D) one. The 2D problem can be obtained by representing
HBARs as planar 2D structure, for example. Here, we provide a brief
description of the HK method. We encourage the interested reader to find a
detailed description of the method elsewhere [65, 66]. The Hankel transform
approach has been mostly used in the field of optics, e.g., Fabry-Perot
cavities [65], however, it has not been applied for acoustics problems, to the
best of our knowledge.

We cast the systems described in Fig. 1 into an axi-symmetric setting about
the $z$-axis. We transform the $x-y$ coordinates into the radial coordinate
$r$, with limits $0\leq r\leq W/2$ discretized into $N$ points as

$r_{j}=W/2\left(\frac{\zeta_{j}}{\zeta_{N}}\right)$ (9)

where $\zeta_{j}$ ($j=1,2,3,..,N$) are the roots of the $J_{1}$ Bessel
function. The round trip operator for the HK approach
$\textbf{RT}_{\text{hk}}$ is given as:

$\textbf{RT}_{\text{hk}}=\textbf{R}_{0}\textbf{P}\textbf{R}_{t_{\text{HBAR}}}\textbf{P}$
(10)

where P, is the $N\times N$ propagator matrix defined as:

P $\displaystyle=(H^{+})^{-1}\tilde{G}H^{+},$ (11a) $\displaystyle H^{+}_{ij}$
$\displaystyle=\frac{W^{2}}{2\zeta_{N}^{2}}\frac{J_{0}(\zeta_{j}\zeta_{j}/\zeta_{N})}{\zeta_{N}^{2}J^{2}_{0}(\zeta_{j})},$
(11b) $\displaystyle\tilde{G}_{ij}$
$\displaystyle=\text{exp}\left(-i\frac{2\zeta^{2}_{j}}{W^{2}}\right)\delta_{jj}.$
(11c)

Here, $H^{+}$ is the Hankel transform whose matrix inverse is taken to be an
approximation of the inverse Hankel transform, and $\tilde{G}$ is the Green’s
function obtained from the Fourier transform of the Fresnel propagator in the
paraxial approximation. $\textbf{R}_{t_{\text{HBAR}}}$, and $\textbf{R}_{0}$
are $N\times N$ diagonal reflection matrix operators at $z=t_{\text{HBAR}}$
and $z=0$ matrix operators, respectively and are defined as

$\displaystyle\textbf{R}_{0,jj}$ $\displaystyle=R_{0,2}(r_{j}),$ (12a)
$\displaystyle\textbf{R}_{t_{\text{HBAR}},jj}$
$\displaystyle=R_{t_{\text{HBAR}},2}(r_{j}).$ (12b)

Since the dimensionality of the problem is reduced from 3D to 2D, the Hankel
eigenvalue problem $\textbf{RT}_{\text{hk}}u(r)=\Lambda u(r)$ can be solved to
obtain the axisymmetric scalar phonon modeshaps $u(r)$ and their corresponding
eigenvalue $\Lambda$. This method, unlike the FBPM/a-FBPM approaches, does not
require a Fox-Li-like iterative process and can also predict higher-order
phonon modes. Although this approach has a limited applicability due to its
axi-symmetric constraints, it was helpful in obtaining phonon modes and
lifetimes of various HBAR structures considered in this study in a
computationally inexpensive manner. The expedited analysis allowed us to
analyze a larger numbers of HBAR structures, as we show in the Results and
Discussion section.

####  II.3.4 Diffraction-Limited Phonon Lifetime

In this article, we estimate the diffraction-limited lifetime of phonons using
three different methods: (A) Eigenvalue method, (B) Exponential curve fitting
method and (C) Clipping method.

(A) Eigenvalue method: In this method, we obtain the eigenvalues using the
adaptive FBPM simulation and use them to estimate the phonon lifetime. The
elastic energy contained in a acoustic beam is given by
$E\propto\int|\textbf{u}|^{2}dxdy$. We do not consider the $z$-dependence of u
for the energy calculation, since the mode shape largely remains unaltered
throughout the thickness. Starting with a cavity phonon mode with a total
elastic energy $E_{tot}$, the ratio of the elastic energy left in the cavity
after one round trip ($E_{in}$) to $E_{tot}$, is given as
$E_{in}/E_{tot}=\Lambda^{2}$. Here, $\Lambda$ is the real eigenvalue obtained
using the a-FBPM simulations. We assume an exponential decay of energy with
propagation time, $E_{in}=E_{tot}e^{-t_{\text{RT}}/\tau}$. Here,
$t_{\text{RT}}=2t_{\text{HBAR}}/v_{2}$ is the time taken by shear waves to
complete a round trip. $\tau$ is the phonon lifetime of the phonon mode and
can then be obtained from

$\tau=\frac{-2t_{\text{HBAR}}}{v_{2}\text{ln}(\Lambda^{2})}.$ (13)

Since these computations are performed on a finite $x-y$ mesh grid, $\Lambda$,
and therefore $\tau$ can be sensitive to the mesh density. Hence, it is
important to select an optimal mesh density to ensure the convergence of
$\Lambda$ since it directly influences $\tau$. In Fig. 5, we show the
variation of $\Lambda$ with mesh density $N_{x}$ ($=N_{y}$) for both the low-
diffraction ($R_{\text{YIG}}=200\mu$m) and the high-diffraction
($R_{\text{YIG}}=10\mu$m) cases. We find that the variation of $\Lambda$ is
much less than 1% for both cases when $N_{x}\geq 1024$ and above.
Consequentially, we choose $N_{x}=N_{y}=1024$ to compute phonon modes and
lifetimes for all cases presented in this article.

Figure 5: Eigenvalues of two planar HBAR structures with
$R_{\text{YIG}}=200\mu$m (red) and $10\mu$m (blue): Variation of simulated
$\Lambda$ for meshes with density ranging from $128\times 128$ to $2400\times
2400$ points, keeping the width $W$ fixed at 1200$\mu$m. The circles with
arrows point to the axis corresponding to each plot.

(B) Exponential curve fitting method: In this method, we allow the initial
beam to travel for multiple round trips and estimate the diffraction loss by
evaluating the overlap of the mode profile with the initial input profile,
after each round trip. We estimate the overlap using the following function:

$I(t)=\frac{|\left<\textbf{u}_{0}(t)|\textbf{u}_{0}(0)\right>|^{2}}{|\left<\textbf{u}_{0}(0)|\textbf{u}_{0}(0)\right>|^{2}}.$
(14)

Here, $t$ is the integral multiple of the time taken for each round trip,
$t_{\text{RT}}$. The propagation is continued until $I(t)<0.1$ is reached or
the time is $t>3$ ms, whichever is reached first. In Fig. 6, we show the decay
of the overlap function, $I(t)$, for beams traveling in a HBAR structure with
$R_{\text{YIG}}=200\mu$m. We show the different decays in this figure to
highlight the effect of the finite width of the transducer (YIG film) on the
propagating beam. The transducer induces a localizing effect on the beam that
can be modeled by introducing a phase factor. We find that $I(t)$ decays at a
much slower rate when the YIG-induced phase factor is considered, compared to
the case without phase factor. We fit the two decays with exponential
functions and obtain the phonon lifetimes as $\tau^{wophase}=0.1$ ms and
$\tau^{wphase}=14.1$ ms, respectively. The remarkable two-orders of magnitude
difference between these two estimates highlights the importance of including
the actuator phase effects for such an analysis. Subsequently, we included the
phase effects in all our analyses to obtain the phonon lifetimes. A past study
has used a similar approach to estimate the lifetime of longitudinal HBAR
phonons in an AlN/Sapphire structure [57]. However, the phonons were treated
as superpositions of Bessel functions in a large simulation window without
accounting for the effect of the localization induced by the AlN transducer.
Interestingly, we find that the lifetime predicted include the phase effects
closely matches with the predicted $\tau$ of $15.5$ ms using the eigenvalue
method, discussed earlier.

Figure 6: Decay of the overlap function, $I(t)$, for propagating beams in a
HBAR structure with $R_{\text{YIG}}=200\mu$m. Solid lines indicate $I(t)$
while dashed lines indicate their respective exponential fits. The finite
width of the YIG film induces a localizing effect on the beam, modeled with a
phase factor. $I(t)$ decays at a much slower rate when the phase effect is
included (blue, top and right axes) compared to the no-phase-effect case (red,
bottom and left axes).

(C) Clipping method: In this method, we obtain
$E_{in}\propto\int_{V_{in}}|\textbf{u}|^{2}dxdy$, contained in a volume,
$V_{in}$, of a cylinder spanning YIG disc and the thickness of the HBAR
structure. We calculate $E_{in}$ for a converged modeshape, obtained with
a-FBPM or Hankel method. With the $E_{in}$ calculated, $\tau$ is obtained as
$\tau=\frac{-2t_{\text{HBAR}}}{v_{2}\text{ln}(E_{in}/E_{tot})}$. This could be
taken as a lower limit to the phonon lifetime since it assumes that all energy
outside the cylinder of interest is lost after a round trip. This method is
similar to the clipping method used in Fabry-Perot cavities [65, 66] and HBAR
resonators excited opto-mechanically using photon beams [73]. In these optical
systems, the clipping method is more applicable since the input photon beam
spans the entire XY plane and is clipped by the localizing finite cross-
section mirrors or domes. In our system, the input acoustic beam is already
laterally confined so this method may have limited applicability. We still
include the discussion here as a consistency check since the method results in
the lower bound for the predicted lifetimes.

###  II.4 Magnon-Phonon Coupling

The magnon-phonon coupling strength, $g_{mb}$, is given by

$g_{mb}=\frac{B}{\sqrt{A_{\nu}Q_{\eta}}}\left|\int_{V_{\text{YIG}}}\left[\frac{du_{x}^{*}}{dz}m_{x}+\frac{du_{y}^{*}}{dz}m_{y}\right]d\textbf{r}\right|.$
(15)

Here, $B=7\times 10^{5}$ J/m3 is the magnetoelastic constant of YIG,
$(m_{x},m_{y})$ are the x-y components of the magnetization
$\textbf{m}=\textbf{m}_{\text{YIG}}$, $(u_{x}^{*},u_{y}^{*})$ and
$(\frac{du_{x}^{*}}{dz},\frac{du_{y}^{*}}{dz})$ represent the complex
conjugates of the x-y components of displacement u and their corresponding
shear strains, respectively. We integrate the coupling between the magnon and
the phonon modes in the volume of the YIG film, $V_{\text{YIG}}$, since the
magnon modes reside in YIG. The term $\sqrt{A_{\nu}Q_{\eta}}$ in the
denominator is a normalizing factor, as we describe below. We normalize the
magnon modes as

$\frac{M_{s}}{\gamma}\int_{V_{\text{YIG}}}\textbf{m}_{\nu}^{*}(\textbf{r})\cdot(\textbf{k}\times\textbf{m}_{\nu^{\prime}}(\textbf{r}))d\textbf{r}=-iA_{\nu}\delta_{\nu,\nu^{\prime}},$
(16)

and the HBAR phonon modes as

$2\omega_{\eta}\rho\int_{V_{\text{HBAR}}}\textbf{u}_{\eta}^{*}(\textbf{r})\cdot\textbf{u}_{\eta^{\prime}}(\textbf{r})d\textbf{r}=Q_{\eta}\delta_{\eta,\eta^{\prime}}.$
(17)

where, $\mu_{0}M_{s}=0.175$T is the saturation magnetization,
$\gamma=2\pi\times 28.5$ GHz/T is the gyromagnetic ratio, $\rho=7080$ kg/m3 is
the GGG material density, k is the unit vector along the z axis, $\nu$ and
$\eta$ are magnon and phonon mode indices, respectively, with $\omega_{\eta}$
being the corresponding phonon mode frequency. Here, $V_{\text{HBAR}}$ is the
volume of the entire HBAR structure. $A_{\nu}$ and $Q_{\eta}$ are phonon and
magnon normalization constants, respectively and are of the same
dimensionality. If zero diffraction is assumed, the acoustic energy is
completely localized in the volume $V_{in}$, resulting in a complete overlap
of lateral mode profiles of magnon and phonon modes. We obtain the zero
diffraction magnon-phonon coupling strength to be $g^{0}_{mb}/{2\pi}=1.13$
MHz, independent of $R_{\text{YIG}}$. This prediction compares well with a
previous experimental result of $1$ MHz [53], however, is slightly higher than
another experimental result of $0.75$ MHz [54]. Although in HBAR structures
some acoustic energy can lie outside $V_{in}$ due to diffraction which can
impact both $\tau$ $g_{mb}$, which we will discuss in the Results and
Discussion section.

##  III Results and Discussion

Here, we discuss the phonon lifetimes and the magnon-phonon coupling in the
HBAR structures shown in Fig. 1, in the diffraction-limited regime. We
investigate the diffraction-limited performance because we are interested in
improving the scalability of hybrid magnonic devices. The diffraction effects
play an increasingly important role as we reduce the radius of the YIG disc as
the reduced aperture increases the diffraction of the acoustic waves
propagating into the GGG region.

Figure 7: Diffraction-limited performance of planar HBAR structures: (a)
Phonon lifetimes, $\tau$, calculated from the eigenvalue (red), exponential
fitting (blue, purple), and clipping methods (green). (b) Ratio between
magnon-phonon coupling in HBAR structures in various diffraction regimes and
that in the zero-diffraction limit, $g_{mb}/g^{0}_{mb}$ (red). Decrease of
$g_{mb}$ is connected to the spread of modeshape in the high-diffraction
regime, as we decrease $R_{\text{YIG}}$. Increasing amount of acoustic energy
spreads into the GGG region, causing reduced overlap between phonon and magnon
mode profiles in the YIG region. The ratio of acoustic energy spreading out to
the total energy, $E_{out}/E_{tot}$ (black), increases in the high-diffraction
regime. (c) Magnon-phonon cooperativity $C$. $R_{\text{YIG}}$ is varied
keeping all other geometric factors fixed. Solid lines correspond to a-FBPM
predictions, while the circles of the same color correspond to the respective
HK predictions.

####  III.0.1 Diffraction-Limited Phonon Lifetime

In Fig. 7 (a), we show the variation of phonon lifetimes in planar HBAR
structures with decreasing YIG radius. We calculate the lifetimes using the
three methods discussed in the Methods section: (A) eigenvalue (red), (B)
exponential fitting (blue, purple), and (C) clipping methods (green). As can
be noted from the Fig. 7, all methods predict that the lifetimes decrease with
decreasing $R_{\text{YIG}}$ due to increased diffraction. We obtain the
highest and the lowest lifetimes to be $\tau=15.5$ms ($Q=\omega\tau=9.57\times
10^{8}$) and $\tau=10.5\mu$s ($Q=6.48\times 10^{5}$), for HBAR structures with
$R_{\text{YIG}}=200\mu$m and $10\mu$m, respectively, using the eigenvalue
method. These results indicate that the lifetimes are reduced by more than
three orders of magnitude in the high-diffraction regime. The eigenvalue
method results in the highest lifetime estimates, among the three methods
considered (red). For the low-to-medium diffraction cases, the lifetimes
predicted by the eigenvalue method (red) match closely with those predicted
from the exponential fitting method (blue), when we include the phase effects
due to the YIG film. This can be explained using the following arguments. The
input beam can be described as a superposition of the eigenmodes of the HBAR
cavity. Once the initial input beam undergoes multiple round trips, only the
fundamental mode is likely to survive, while the rest of the mode components
decay at a faster rate. When we operate at the overtone frequencies of one of
the fundamental modes, the fundamental mode becomes the dominant mode. It can
be argued that this dominant mode is the same as the fundamental eigenmode
found by solving the eigenvalue equation, Eq. 8. Therefore, the decay of the
overlap function of the input field is expected to have a similar character to
that of the decay of the dominant eigenmode that it is composed of. As a
result, we obtain similar lifetimes values using the two different methods.

However, we observe 35-fold reduced lifetimes predictions from the exponential
fit compared to the eigenvalue method, for the highest diffraction case
($R_{\text{YIG}}=10\mu$m). In the exponential fit method, we perform the
exponential fitting of the acoustic overlap function, $I(t)$, starting at
$t=0$. Due to high diffraction, significant amount of the energy of the input
acoustic beam spreads out laterally even before the YIG film introduces
localizing effects on the propagating beam. This results in a rapid reduction
of $I(t)$ during the first few round trips. The loss of acoustic overlap
during the initial transient phase results in lower $\tau$ predictions. We
expect that the two lifetime predictions will be closer if the exponential fit
is obtained after a steady state has been achieved. Note that these
exponential fitting predictions are computed by including the phase effects
due to the YIG film. We also show the lifetimes predicted from the exponential
fitting method, when the phase effects due to the YIG film are ignored
(purple). These effects were often ignored in past analyses of HBAR phonons,
however, as we show in Fig. 7, they significantly affect the predictions. We
obtain the $\tau^{wophase}$’s that are consistently lower than the
$\tau^{wphase}$’s, as also discussed in the Methods section (Fig. 6). The
difference between the two predictions $\tau^{wophase}$ and $\tau^{wphase}$
decreases monotonously from 135 to 1.35 fold as we decrease $R_{\text{YIG}}$
from $=200\mu$m to $=10\mu$m, respectively. This is expected since the effect
of the YIG film width on the mode diminishes as we approach
$R_{\text{YIG}}\rightarrow 0$.

We also show the lifetime predictions from the clipping method (green) which
assumes that the energy outside the cylindrical volume of interest, $V_{in}$,
is completely lost after a round trip. We find the lifetimes to be more than
an order of magnitude lower than those predicted from the eigenvalue method
for all cases considered. We find the lifetimes are an order of magnitude
lower than those predicted from the eigenvalue method for all cases
considered. We obtain the lifetime in the highest diffraction regime to be
$\tau=0.311\mu$s and the respective quality factor, $Q$ is $1.92\times
10^{4}$. Interestingly, these values are still marginally greater than
experimentally observed values of 0.25$\mu$s and $1.45\times 10^{4}$,
respectively [54], for a device operating in the low-diffracting regime. The
reason for this discrepancy is that the experimental values are obtained for
HBAR devices operating at room temperatures. In this limit, the performance of
HBAR structures is limited by the material attenuation effects [55, 56]. In
our study, we ignore the material attenuation effects and only discuss the
diffraction-limited performance. Our predictions will be more representative
for devices operating in the cryogenic or potentially milli Kelvin regimes.
For example, the HBAR devices could be envisioned to effectively couple with
superconducting qubit systems, etc., which typically operate in the milli
Kelvin regime.

We like to address here the uncertainties associated with the numerical
prediction of the HBAR phonon lifetimes using different methods. In the
eigenvalue method, we obtain the lifetimes from the ratio between the HBAR
thickness ($t_{\text{HBAR}}$) and a function of the phonon velocity and the
eigenvalue, $\Lambda$, as shown in Eq. 13. Following Eq. 13, an error
propagation relation can be written as

$\left|\frac{\text{d}\tau}{\tau}\right|=\left|\frac{v_{2}\tau}{t_{\text{HBAR}}}\times\frac{\text{d}\Lambda}{\Lambda}\right|\propto\left|\tau\times\frac{\text{d}\Lambda}{\Lambda}\right|.$
(18)

Here, the order of magnitude of the prefactor
$\frac{v_{2}}{t_{\text{HBAR}}}\sim 10^{7}$. However, the prefactor does not
affect the predictions in our study since we keep the material and thickness
fixed. Equation 18 suggests that the uncertainty
$\left|\frac{\text{d}\tau}{\tau}\right|$ increases monotonously with $\tau$
when other factors are fixed. It can be deduced from Equation 18 that to
predict a lifetime within $x\%$ of variability, $\Lambda$ has to be predicted
to be within $\frac{t_{\text{HBAR}}x}{v_{m}\tau}\%$ of variability. This
implies that we have to be increasingly more stringent with our $\Lambda$
convergence criteria as $\tau$ increases. Figure 8 shows that the conditions
for desired accuracy (relative tolerance) needed for $\Lambda$ predictions
become more stringent for high-$\tau$ (high-$Q$) systems. When $\tau=0.1\mu$s
and the desired accuracy is set to $10\%$, $\Lambda$ must be predicted within
$15\%$ accuracy, which is attainable following our numerical procedure. The
accuracy requirement increases fast when the lifetime values are much higher.
For example, when $\tau=1$ms, and the desired accuracy is set to $10\%$ of
$\tau$ value, $\Lambda$ must be predicted within ${\sim}10^{-3}\%$ of
accuracy. For our HBAR with $R_{\text{YIG}}=200\mu$m system, we obtain
$\left|\frac{\text{d}\tau}{\tau}\right|=7.9\times 10^{-4}$ with
$\left|\frac{\text{d}\Lambda}{\Lambda}\right|=4.7\times 10^{-7}$ between last
two restarts. On the other hand for the HBAR with $R_{\text{YIG}}=10\mu$m
system, $\left|\frac{\text{d}\tau}{\tau}\right|=3.9\times 10^{-7}$ with
$\left|\frac{\text{d}\Lambda}{\Lambda}\right|=3.8\times 10^{-8}$. However,
such high accuracy is challenging to achieve numerically because of the
computational expense involved in simulating a large number of iterations. We
discuss an alternative HK approach to expedite the numerical analysis.

Figure 8: Relative tolerance of $\Lambda$ at various $\tau$ required to
predict $\tau$ with $10\%$ accuracy.

####  III.0.2 Magnon-Phonon Coupling in Diffraction-Limited Regime

We now turn to discuss the effect of diffraction on the magnon-phonon coupling
strength, $g_{mb}$. We compute the ratio between coupling strength, $g_{mb}$,
in the planar HBAR structures, and that in the zero-diffraction limit,
$g^{0}_{mb}$, to estimate the effect of diffraction on magnon-phonon coupling.
In Fig. 7(b), we show the variation of the ratio between coupling strengths,
$g_{mb}/g^{0}_{mb}$, in planar HBAR structures with decreasing YIG radius. The
ratio is equal to 1 for low-diffraction cases with $R_{\text{YIG}}\gtrsim
100\mu$m, as expected. As we reduce $R_{\text{YIG}}$ to $40\mu$m, $g_{mb}$
only changes by 5% from the $g^{0}_{mb}$ limit. Note that this case typically
corresponds to a strongly diffracting regime, however, the effect of
diffraction is minimized. This is because the YIG disc helps to localize the
beam and thus, preserve the magnon-phonon overlap and the coupling strength.
However, we observe a sharp decline of $g_{mb}$ for $R_{\text{YIG}}\leq
40\mu$m. $g_{mb}$ drops to $77.4\%$ of $g^{0}_{mb}$ for HBAR with
$R_{\text{YIG}}=20\mu$m and $55.7\%$ of $g^{0}_{mb}$ for
$R_{\text{YIG}}=10\mu$m, the highest diffraction case considered here. The
decrease of $g_{mb}$ can be explained in the following way. In the zero
diffraction limit, the acoustic energy of a propagating beam is completely
localized in the volume $V_{in}$, of a cylinder with radius equal to
$R_{\text{YIG}}$ and the thickness of the HBAR structure. This results in a
complete overlap of lateral mode profiles of magnon and phonon modes. However,
as we decrease $R_{\text{YIG}}$, the modeshape spreads out beyond $V_{in}$,
due to diffraction. Correspondingly, some acoustic energy also leaks out of
$V_{in}$ and spreads into the GGG region. As a result, there is reduced
overlap between the phonon and the magnon mode profiles in the YIG region. We
refer to the energy of the mode lying outside $V_{in}$ as $E_{out}$. We show
the change of the ratio $E_{out}/E_{tot}$ with decreasing $R_{\text{YIG}}$ in
Fig. 7(b) (black). As we can see the ratio of acoustic energy of the eigen
mode residing outside $V_{in}$ to the total energy increases in the high-
diffraction regime.

Finally, we combine the phonon lifetimes and the magnon-phonon coupling
strength to determine the performance figure of merit of the HBAR structures,
defined by the magnon-phonon cooperativity,
$C=4g^{2}_{mb}/\kappa_{m}\kappa_{b}=4g^{2}_{mb}\tau_{m}\tau_{b}$. We show the
variation of $C$ with $R_{\text{YIG}}$ in Fig. 7 (c). Using $\tau$ predicted
from the eigenvalue method, we obtain a monotonically decreasing diffraction-
limited $C$ ranging from $21.9\times 10^{4}$ to $46.4$, as we decrease
$R_{\text{YIG}}$ from $200\mu$m to $10\mu$m. On the other hand, using $\tau$
predicted from the clipping method, $C$ is obtained in the range between
$1.2\times 10^{4}$ and $1.4$. Note that even the HBAR with
$R_{\text{YIG}}=10\mu$m has a high cooperativity. This implies that if
material acoustic attenuation effects are eliminated (e.g. operating in the
milli Kelvin regime), one can achieve remarkable scalability with planar HBAR
structures. However, the performance of the HBAR structure with
$R_{\text{YIG}}=10\mu$m is limited by a lifetime of $10.5\mu$s ($0.31\mu$s)
predicted by the eigenvalue (clipping) methods and only $55.7\%$ of maximum
coupling strength $g^{0}_{mb}$. Thus, it is desirable to pursue design
approaches to further improve these performance parameters.

In addition to the a-FBPM predictions discussed so far, we also show
predictions from HK method in Fig. 7 using circles, with colors matching to
their corresponding a-FBPM predictions. Note that we use $|\Lambda|$ instead
of $\Lambda$ in Eq. 13 to calculate $\tau$ using the HK method. As can be
noted from Fig. 7, $\tau$, $g_{mb}$, and $C$ predicted by the HK method are in
an excellent agreement with a-FBPM predictions, for all $R_{\text{YIG}}$ cases
considered. The close match validates our implementation of the HK method. As
we have discussed earlier, a-FBPM has broader applicability compared to HK
method, and can be applied for systems with anisotropy in material properties
and transducer shapes. However, we find that it is an expensive and uncertain
process to converge the residual R as the rates of convergence for different
HBAR structures are slow and variable. Also, one needs to select the number of
round trips before restarting the simulation using a trial-and-error approach.
In comparison, HK method is computationally inexpensive to make predictions.
Additionally, it can be applied to the YIG/GGG HBAR systems since they are
isotropic with YIG transducer being a circular thin film disc. Henceforth, we
only present predictions from the HK method in this article.

####  III.0.3 Enhancing Scalability and Performance in Diffraction-Limited
Regime

The discussions presented in the previous subsections show that the phonon
lifetimes are reduced by three orders or magnitude in the high-diffraction
regime, when $R_{\text{YIG}}=10\mu$m, compared to low-diffraction cases. The
magnon-phonon coupling is also reduced by a factor of two. The phonon lifetime
reduction plays a dominant role and results in a three orders of magnitude
reduction of the cooperativity, the performance figure of merit of the HBAR
structures. Here, we discuss a strategy to mitigate the diffraction losses and
improve the performance of the HBAR structures. We illustrate that the use of
focusing dome-like surface structures could significantly improve the
performance figure of merit of HBAR structures. Confocal HBAR structures have
been demonstrated to have Q-factors on the order of $10^{7}$ while resulting
in $10^{3}$-fold reduction in device volumes [73]. However, they have not been
employed in hybrid magnomechanical systems. Our results show that a confocal
geometry could enhance the phonon lifetimes and improve $g_{mb}$ due to
lowered mode volume of phonons resulting from reduced lateral spread of the
phonon modes.

We first discuss a strategy to identify the shape of the dome that will result
in the highest $g_{mb}$. We find that $g_{mb}$ is highest when there is
maximum overlap between the phonon and the magnon modeshapes. We assume that
the CHBAR phonons can be described by Laguerre-Gaussian (LG) functions. The
fundamental mode, LG00, is of the form: $u_{y}\propto e^{-R^{2}/w_{0}^{2}}$
due to the axi-symmetry of our system. Here, $R$ is the radial coordinate
$\sqrt{x^{2}+y^{2}}$. $w_{0}$ is the waist of the Gaussian beam at the $z=0$
of the HBAR structure and is an adjustable parameter that can be varied by
changing the dome radius of curvature $R_{\text{curv}}$, $t_{\text{HBAR}}$, or
$\omega$. We ignore the effects of lateral extent of the structure and the
dome, and the localizing effects of the YIG film. We obtain theoretical
estimates of the overlap of the CHBAR phonon modeshapes with Bessel magnon
modes and extract the optimal dome shape using the condition where the overlap
is maximum. We normalize the modeshapes according to Eqs. 16, 17 and calculate
the coupling strengths using Eq. 15. We show our predicted values of the
$g_{mb}/g^{0}_{mb}$ ratio in Fig. 9, as a function of the beam waist, $w_{0}$.
We find that $g_{mb}/g^{0}_{mb}$ ratio peaks at $w_{0}/R_{\text{YIG}}~{}\sim
0.65$. Interestingly, we find that this optimal ratio, $w_{0}/R_{\text{YIG}}$,
is constant and does not depend on $R_{\text{YIG}}$, for all the CHBAR
structures considered in this article. If we decrease $w_{0}$ lower than
$w_{0}~{}\sim 0.65R_{\text{YIG}}$, the coupling strength decline sharply. On
the other hand, the coupling decrease slowly if $w_{0}$ is increased beyond
the optimal value. We use this optimal $w_{0}$ value to provide an initial
estimate of the radius of curvature, $R_{\text{curv}}$, of the dome [74]:

$\displaystyle R_{\text{curv}}$
$\displaystyle=\frac{1}{Re[1/(q_{0}+t_{\text{HBAR}})]},$ (19a)
$\displaystyle\text{with }q_{0}$ $\displaystyle=\frac{i\pi
w_{0}^{2}}{\lambda}.$ (19b)

This expression indicates that we can get multiple estimates of
$R_{\text{curv}}$ by varying $t_{\text{HBAR}}$, $\lambda$ and $w_{0}$.

Figure 9: Determination of dome-shape that results in the highest magnon-
phonon coupling, $g_{mb}$, in CHBAR structures: Variation of
$g_{mb}/g^{0}_{mb}$ with the waist of the fundamental phonon mode.
$g_{mb}/g^{0}_{mb}$ is maximum when $w_{0}/R_{\text{YIG}}~{}\sim 0.65$.
Figure 10: Performance of CHBAR structures with varied radius of curvature,
$R_{\text{curv}}$, of the dome-shape: (a) Phonon lifetimes, $\tau$, calculated
from the eigenvalue (red) and clipping methods (green), (b) Ratio between
magnon-phonon coupling in HBAR structures in various diffraction regimes and
that in the zero-diffraction limit, $g_{mb}/g^{0}_{mb}$ and (c) Magnon-phonon
cooperativity, $C$, calculated using the HK method. $R_{\text{curv}}$ is
varied keeping $R_{\text{YIG}}=10\mu$m and $R_{\text{cross}}=60\mu$m fixed.
Solid lines represent corresponding values for a planar HBAR structure with
$R_{\text{YIG}}=10\mu$m, while the circles of the same color correspond to the
CHBAR values.

Using Eqs. 19 (a-b), we obtain the analytical estimate to be
$R_{\text{curv}}/t_{\text{HBAR}}{\sim}1.5$ that corresponds to
$w_{0}/R_{\text{YIG}}~{}{\sim}0.65$ and the peak coupling value as shown in
Fig. 9. Note that this $R_{\text{curv}}$ estimate is obtained only to maximize
$g_{mb}/g^{0}_{mb}$, and the analysis does not consider the effects of the
dome shape on phonon lifetimes and cooperativities. Also, this analysis
ignores YIG localization effects and assumes phonon modeshapes to be perfectly
Gaussian. To obtain a better picture of the overall performance of these CHBAR
structures, we carry out a numerical study using the HK method. We include the
effects of the lateral extents of the dome and the device, the attenuation
window, and the localization effects of the YIG film, which are not considered
for the analytical predictions. We modify the equations for analyzing CHBAR
structures in the following way. We modify Eq. 6 and the radial form of Eq.
11c to include the phase induced by the dome surface as

$R_{t_{\text{CHBAR}},m}(x,y)=\begin{cases}e^{ik_{z0,m}r^{2}/R_{\text{curv}}},&\text{if
}r\leq R_{\text{cross}}\\\ 1,&\text{if
}R_{\text{cross}}<r\leq\frac{W_{\text{eff}}}{2}\\\
0.&\text{otherwise}\end{cases}$ (20)

We sweep across various values centered around the analytical estimate,
$R_{\text{curv}}/t_{\text{HBAR}}{\sim}1.5$, to identify the optimal
$R_{\text{curv}}$ for the CHBAR structures. We choose the CHBAR structure with
$R_{\text{YIG}}=10\mu$m and the dome $x-y$ cross-section radius,
$R_{\text{cross}}=60\mu$m, to perform the analysis.

We show the phonon lifetime, magnon-phonon coupling and cooperativity in Fig.
10, to highlight the effect of the focusing dome on the performance of CHBAR
structure with $R_{\text{YIG}}=10\mu$m. We show the variation of $\tau$ with
$R_{\text{curv}}/t_{\text{HBAR}}$ in Fig. 10(a). The predicted $\tau$ shows a
peak at $R_{\text{curv}}/t_{\text{HBAR}}=1.09$, a value less than the
analytical prediction for maximum $g_{mb}$. The peak value reaches
$\tau{\sim}218$ms using the eigenvalue method (red circles). Incidentally,
this value is ${\sim}$500 times larger than the value at
$R_{\text{curv}}/t_{\text{HBAR}}{\sim}1.5$, justifying the necessity of
performing the sweep. The $\tau$ peak value for the CHBAR structure obtained
using the eigenvalue (clipping) method is ${\sim}1.7\times 10^{4}$
(${\sim}7.3\times 10^{4}$) factor greater than the corresponding planar HBAR
$\tau$ prediction, shown with red horizontal line in Fig. 10(a). We show the
predictions from the clipping method using green circles. To obtain these
predictions, we take $V_{in}$ to be the volume defined by the dome’s lateral
cross-sectional area and the structure thickness. The predicted $\tau$ peaks
at ${\sim}22.7$ms, lower than those predicted by the eigenvalue method, as
expected. We note a sharp reduction in $\tau$ is observed for
$R_{\text{curv}}/t_{\text{HBAR}}<1$. This corresponds to the situation when
the center of curvature of the dome is inside the thickness of the HBAR
structure, resulting in a negative ‘$g-$parameter’, given by
$g=1-\frac{t_{\text{HBAR}}}{R_{\text{curv}}}$. It has been shown that one
cannot obtain real and finite Gaussian beam solutions for a negative
$g-$parameter [73]. We find that this is the case for the modes corresponding
to $R_{\text{curv}}/t_{\text{HBAR}}<1$, that is, they are lossy modes which do
not have well-defined Gaussian characters.

We show the corresponding $g_{mb}/g^{0}_{mb}$ estimates in Fig. 10(b). We
obtain a peak value of 0.99 at $R_{\text{curv}}/t_{\text{HBAR}}{\sim}1.4$, a
value close to the analytical estimate of ${\sim}1.5$. However, $\tau$ is
reduced by more than two orders of magnitude at this point. A key aspect to
note for these systems is that reducing the phonon mode volume by means of
focusing does not necessarily improve $g_{mb}$ beyond its maximum achievable
value, $g^{0}_{mb}$, which represents the case when there is complete lateral
overlap of phonon and magnon modes. We use the $\tau$ and $g_{mb}$ values to
calculate the figure of merit, $C$, for all the CHBAR structures considered.
As can be seen from Fig. 10(c), the $C$ predictions reach peak value at
$R_{\text{curv}}/t_{\text{HBAR}}=1.09$. This further illustrates that $\tau$
dominates the performance of these structures, since there is a ceiling to the
achievable $g_{mb}$. The peak $C$ is calculated to be $2.5\times 10^{6}$
($2.6\times 10^{5}$) using the eigenvalue (clipping) method. This CHBAR
structure shows several orders of magnitude improvements in $\tau$ and $C$
compared to its planar counterpart, and reaches up to $90\%$ of $g^{0}_{mb}$
with a scalability factor of $S~{}{\sim}64$, limited by the $x-y$ cross-
sectional area of the dome. Note that we only vary $R_{\text{curv}}$ for this
analysis keeping the $R_{\text{cross}}$ fixed at $60\mu$m which limits the
lateral scalability of these devices.

Figure 11: Performance of CHBAR structures with varied radius of cross-
section, $R_{\text{cross}}$, of the dome-shape: (a) Phonon lifetimes, $\tau$,
calculated from the eigenvalue (red) and clipping methods (green), (b) Ratio
between magnon-phonon coupling in HBAR structures in various diffraction
regimes and that in the zero-diffraction limit, $g_{mb}/g^{0}_{mb}$ and
(c)Magnon-phonon cooperativity, $C$, calculated using the HK method.
$R_{\text{cross}}$ is varied keeping $R_{\text{YIG}}=10\mu$m and
$R_{\text{curv}}/t_{\text{HBAR}}=1.09$ fixed. Solid lines represent
corresponding values for a planar HBAR structure with $R_{\text{YIG}}=10\mu$m,
while the circles of the same color correspond to the CHBAR values.

We perform additional analysis to further improve the scalability of the CHBAR
devices. We reduce the cross-sectional area of the dome while keeping
$R_{\text{curv}}/t_{\text{HBAR}}=1.09$ and $R_{\text{YIG}}=10\mu$m fixed. We
show the variation of $\tau$, $g_{mb}$, and $C$ with
$R_{\text{cross}}/R_{\text{YIG}}$ in Fig. 11. As we decrease
$R_{\text{cross}}/R_{\text{YIG}}$ from $6$ to $1$, the phonon lifetime is
reduced by more than 4 orders of magnitude from ${\sim}218$ ms (22.7 ms) to
$4.83\mu$s ($0.134\mu$s) as predicted by the eigenvalue (clipping) method. The
magnon-phonon coupling shows a sharp decline for
$R_{\text{cross}}/R_{\text{YIG}}<3$ reaching the lowest value of
$g_{mb}/g^{0}_{mb}~{}{\sim}0.32$ at $R_{\text{cross}}/R_{\text{YIG}}=1$. While
the coupling strengths are minimally affected for
$R_{\text{cross}}/R_{\text{YIG}}\geq 3$. The cooperativity, $C$, vs
$R_{\text{cross}}/R_{\text{YIG}}$ follows a similar trend to that of $\tau$,
as shown in Fig. 11(c). Using the predictions from the eigenvalue (clipping)
method, we obtain that $C$ is decreased by more than 5 orders of magnitude
from $2.5\times 10^{6}$ ($2.6\times 10^{5}$) to 7.2 (0.19) as we decrease
$R_{\text{cross}}/R_{\text{YIG}}$. Additionally, we find that the performance
of these CHBAR structures are lowered below those of planar HBAR values, shown
with solid lines, if $R_{\text{cross}}/R_{\text{YIG}}<2$. Considering all the
different aspects, we find that the optimal geometry for the CHBAR structures
is represented by the condition $R_{\text{cross}}/R_{\text{YIG}}=4.5$. The
phonon lifetime is long and the lifetime variation is relatively flat in the
range between 4.5 and 6. The numerical analysis of this CHBAR structure
($R_{\text{curv}}/t_{\text{HBAR}}=1.09$, $R_{\text{YIG}}=10\mu$m,
$R_{\text{cross}}/R_{\text{YIG}}=4.5$) using the eigenvalue (clipping) method
predicts that the diffraction-limited lifetime of $\tau=$144.7 ms (11.1 ms)
and a scalability of $S=113$ could be achieved at a cooperativity
$C=1.64\times 10^{6}$ ($1.25\times 10^{5}$). The analysis and the results
presented in Figs. 10 and 11 provide a proof-of-concept that a focusing dome
improves the performance of hybrid magnonic HBAR structures, for quantum
memory and transduction applications.

##  IV Conclusion and Outlook

In aummary, our study provides key insights into the diffraction-limited
performance of YIG/GGG HBAR hybrid magnonic devices. Additionally, we
establish a prediction approach that can be used for designing scalable, long-
lasting quantum memories and efficient quantum transduction systems. We
present analytical and numerical analyses of the diffraction-limited BAW
phonon lifetimes, modeshapes, and the magnon-phonon coupling strengths in
planar and confocal HBAR structures. We use (1) Fourier beam propagation
(FBPM) and (2) Hankel transform (HK) eigenvalue problem approaches to analyze
the shear wave phonon modes in the HBAR structures. The FBPM approach has been
widely used to analyze beam propagation in the field of optics, and more
recently, to study HBAR phonons in planar and confocal HBAR (CHBAR)
structures. Here, we implement a reformulated approach that allows us to
achieve a seven-fold speed up of the computation time. However, we find that
the number of round trips needed in FBPM method to convergence the phonons are
often large and unpredictable. We implement an adaptive FBPM (a-FBPM) approach
that mostly allows to overcome the slow convergence issues over the standard
approach. The a-FBPM still suffers from convergence challenges for the
confocal HBAR structures. We implement the HK method that leverages the axi-
symmetry of the YIG/GGG HBAR system and is computationally inexpensive. The HK
approach has been mostly used in the field of optics, e.g., Fabry-Perot
cavities, however, it has not been applied for acoustics problems, to the best
of our knowledge.

Our analysis predicts that the diffraction-limited $\tau$ of a planar HBAR
structure with lateral YIG dimension, $R_{\text{YIG}}=200\mu$m, is on the
order of milliseconds. A recent study reported that the room temperature
performance of a YIG/GGG HBAR structure is limited by the phonon lifetime at
$0.25\mu$s [54]. The HBAR structure had larger YIG lateral area $0.72$ mm2
than that of the HBARs considered in our study. The difference in $\tau$
values implies that the performance of the previously studied system is not
limited by diffraction losses since diffraction effects are less dominant for
larger YIG lateral area. Instead, the performance is likely to be limited by
material losses. Assuming that the material-limited lifetime could be pushed
to $~{}{\sim}0.1$ ms at mK temperatures, we find that the planar HBAR
structures are not affected by diffraction effects even at
$R_{\text{YIG}}=50\mu$m, which already offers significant 50-fold scalability.
The scalability can be further improved by scaling down the YIG film lateral
area. Increased lifetimes have direct implications on the storage times of
quantum states for quantum memory applications, whereas scalability is
naturally desired to accommodate multiple on-chip memory centers or to have
other on-chip hardware alongside the HBAR structures. We acknowledge that a
full analysis of the material-limited lifetimes is still necessary to obtain a
comprehensive understanding of the shear waves in YIG and GGG material
systems. This may include a clear identification of the operation regimes
(e.g., Akhiezer and Landau-Rumer) at various temperatures and frequencies of
interest.

The performance of planar HBAR devices is heavily dependent on maintaining
perfectly parallel surfaces. We illustrate that the focusing dome structure
allows to mitigate the sensitivity to parallelism and significantly improve
the phonon lifetime and scalability of these systems. To identify the shape of
the dome structure, we theoretically estimate overlap of the CHBAR phonon
modeshapes with magnon modes and extract the shape paramters using the
condition where the overlap is maximum. We refine the estimate with numerical
analysis and obtain a radius of curvature and cross-section radius for which
both $\tau$ and $C$ are optimal and the coupling is close to its peak value.
Overall, we find that ultra-high, diffraction-limited, cooperativities and
phonon lifetimes on the order of ${\sim}10^{5}$ and ${\sim}10$ ms,
respectively, could be achieved using a CHBAR structure with
$R_{\text{YIG}}=10\mu$m. In addition to enhanced $\tau$ and $C$, the confocal
HBAR structure is predicted to provide more than 100-fold improvement of
scalability.

Our study provides key insights into the diffraction-limited performance
induced by the lateral area of the YIG film and the dome structure. Various
other parameters, such as the thickness of the YIG and GGG regions,
geometrical misalignment, and imperfections can also play an important role in
determining the performance. Additionally, the photon and magnon lifetimes and
their coupling strength could also play important role in the overall device
performance. We assume that these parameters remain invariant in our study.
However, we acknowledge that a comprehensive analysis is necessary to
understand to optimize the device performance. Such an analysis could be
performed using the current state-of-the-art machine learning techniques,
which is a promising research direction for the future. Anharmonic phonon
interactions at the surface and interfaces, and the coupling of bulk acoustic
waves to surface acoustic waves, is another interesting direction to explore.
The current study only considers the coupling of the fundamental magnon mode
to the fundamental phonon mode for the given structure. However, it is
possible that the higher-order mode coupling could play a major role
considering the broadband nature of the phonon modes and improved lifetimes
introduced by the dome structure. It will be interesting to explore the
applicability of YIG/GGG HBAR structures for quantum transduction applications
by coupling with other quantum information carriers such as superconducting
qubits, which also operate in the microwave frequency regime and at milli
Kelvin temperatures. Our study is concerned with coherent states, and it will
be important to explore how the insights provided here translate to systems
dealing with non-classical states, e.g. squeezed states, cat states, and Fock
states.

##  V Acknowledgements

We are indebted to Prof. Xufeng Zhang for sharing valuable insights regarding
the outstanding challenges of hybrid magnonic systems and several helpful
discussions. We gratefully acknowledge funding from the Quantum Explorations
in Science & Technology (QuEST) grant provided by the CU Boulder Research &
Innovation Office (RIO) in partnership with the College of Engineering and
Applied Science, the College of Arts and Sciences, JILA, and the National
Institute of Standards and Technology (NIST). We acknowledge the computing
resources provided the RMACC Summit supercomputer, which is supported by the
National Science Foundation (awards ACI-1532235 and ACI-1532236), the
University of Colorado Boulder, and Colorado State University. The Summit
supercomputer is a joint effort of the University of Colorado Boulder and
Colorado State University.

## References

  * Lachance-Quirion  _et al._ [2019] D. Lachance-Quirion, Y. Tabuchi, A. Gloppe, K. Usami, and Y. Nakamura, Hybrid quantum systems based on magnonics, Applied Physics Express 12, 070101 (2019).
  * Li  _et al._ [2020] Y. Li, W. Zhang, V. Tyberkevych, W.-K. Kwok, A. Hoffmann, and V. Novosad, Hybrid magnonics: Physics, circuits, and applications for coherent information processing, Journal of Applied Physics 128 (2020).
  * Huebl  _et al._ [2013] H. Huebl, C. Zollitsch, J. Lotze, F. Hocke, M. Greifenstein, A. Marx, R. Gross, and S. Goennenwein, High cooperativity in coupled microwave resonator ferrimagnetic insulator hybrids, Phys. Rev. Lett. 111, 127003 (2013).
  * Zhang  _et al._ [2014] X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Strongly coupled magnons and cavity microwave photons, Phys. Rev. Lett. 113, 156401 (2014).
  * Tabuchi  _et al._ [2014] 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).
  * Bai  _et al._ [2015] L. Bai, M. Harder, Y. P. Chen, X. Fan, J. Q. Xiao, and C.-M. Hu, Spin pumping in electrodynamically coupled magnon-photon systems, Phys. Rev. Lett. 114, 227201 (2015).
  * Goryachev  _et al._ [2014] M. Goryachev, W. G. Farr, D. L. Creedon, Y. Fan, M. Kostylev, and M. E. Tobar, High-cooperativity cavity qed with magnons at microwave frequencies, Phys. Rev. Appl. 2, 054002 (2014).
  * Zhang  _et al._ [2016a] X. Zhang, N. Zhu, C.-L. Zou, and H. X. Tang, Optomagnonic whispering gallery microresonators, Phys. Rev. Lett. 117, 123605 (2016a).
  * Osada  _et al._ [2016] A. Osada, R. Hisatomi, A. Noguchi, Y. Tabuchi, R. Yamazaki, K. Usami, M. Sadgrove, R. Yalla, and Y. Nomura, M.and Nakamura, Cavity optomagnonics with spin-orbit coupled photons, Phys.Rev.Lett 116, 223601 (2016).
  * Haigh  _et al._ [2016] J. A. Haigh, A. Nunnenkamp, A. J. Ramsay, and A. J. Ferguson, Triple-resonant brillouin light scattering in magneto-optical cavities, Physical Review Letters 117, 133602 (2016).
  * Zhang  _et al._ [2016b] X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Cavity magnomechanics, Science Advances 2, e1501286 (2016b).
  * Hisatomi  _et al._ [2016] R. Hisatomi, A. Osada, Y. Tabuchi, T. Ishikawa, A. Noguchi, R. Yamazaki, K. Usami, and Y. Nakamura, Bidirectional conversion between microwave and light via ferromagnetic magnons, Phys. Rev. B 93, 174427 (2016).
  * Zhu  _et al._ [2020a] N. Zhu, X. Zhang, X. Han, C. Zou, C. Zhong, C. Wang, L. Jiang, and H. Tang, Waveguide cavity optomagnonics for broadband multimode microwave-to-optics conversion, Optica 7, 1291 (2020a).
  * Tabuchi  _et al._ [2015] 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).
  * Lachance-Quirion  _et al._ [2020] 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).
  * Wolski  _et al._ [2020] S. P. Wolski, D. Lachance-Quirion, Y. Tabuchi, S. Kono, A. Noguchi, K. Usami, and Y. Nakamura, Dissipation-based quantum sensing of magnons with a superconducting qubit, Phys.Rev.Lett 125, 117701 (2020).
  * Flower  _et al._ [2019] G. Flower, J. Bourhill, M. Goryachev, and M. E. Tobar, Broadening frequency range of a ferromagnetic axion haloscope with strongly coupled cavity–magnon polaritons, Phys. Dark Universe 25, 100306 (2019).
  * Trickle  _et al._ [2020] T. Trickle, Z. Zhang, and K. M. Zurek, Detecting light dark matter with magnons, Phys. Rev. Lett. 124, 201801 (2020).
  * Wang  _et al._ [2019a] Y.-P. Wang, J. W. Rao, Y. Yang, P.-C. Xu, Y. S. Gui, B. M. Yao, J. Q. You, and C.-M. Hu, Nonreciprocity and unidirectional invisibility in cavity magnonics, Physical Review Letters 123, 127202 (2019a).
  * Zhang  _et al._ [2020a] X. Zhang, A. Galda, X. Han, D. Jin, and V. M. Vinokur, Broadband nonreciprocity enabled by strong coupling of magnons and microwave photons, Phys. Rev. Appl. 13, 044039 (2020a).
  * Zhu  _et al._ [2020b] N. Zhu, X. Han, C.-L. Zou, M. Xu, and H. X. Tang, Magnon-photon strong coupling for tunable microwave circulators, Phys. Rev. A 101, 043842 (2020b).
  * Zhang  _et al._ [2017] D. Zhang, X.-Q. Luo, Y.-P. Wang, T.-F. Li, and J. You, Observation of the exceptional point in cavity magnon-polaritons, Nat.Commun 8, 1368 (2017).
  * Harder  _et al._ [2017] M. Harder, L. Bai, P. Hyde, and C.-M. Hu, Topological properties of a coupled spin-photon system induced by damping, Phys. Rev. B 95, 214411 (2017).
  * Zhang  _et al._ [2019] X. Zhang, K. Ding, X. Zhou, J. Xu, and D. Jin, Experimental observation of an exceptional surface in synthetic dimensions with magnon polaritons, Phys. Rev. Lett. 123, 237202 (2019).
  * Bozhko  _et al._ [2017] D. A. Bozhko, P. Clausen, G. A. Melkov, V. S. L’vov, A. Pomyalov, V. I. Vasyuchka, A. V. Chumak, B. Hillebrands, and A. A. Serga, Bottleneck accumulation of hybrid magnetoelastic bosons, Phys. Rev. Lett. 118, 237201 (2017).
  * Lisenkov  _et al._ [2019] I. Lisenkov, A. Jander, and P. Dhagat, Magnetoelastic parametric instabilities of localized spin waves induced by traveling elastic waves, Phys. Rev. B 99, 184433 (2019).
  * Verba  _et al._ [2019] R. Verba, V. Tiberkevich, and A. Slavin, Wide-band nonreciprocity of surface acoustic waves induced by magnetoelastic coupling with a synthetic antiferromagnet, Phys. Rev. Appl. 12, 054061 (2019).
  * Shah  _et al._ [2020] P. Shah, D. Bas, I. Lisenkov, A. Matyushov, N. Sun, and M. Page, Giant nonreciprocity of surface acoustic waves enabled by the magnetoelastic interaction, Sci Adv 6, eabc5648 (2020).
  * Kittel [1958] C. Kittel, Interaction of spin waves and ultrasonic waves in ferromagnetic crystals, Phys.Rev. 110, 836 (1958).
  * Schlomann [1960] E. Schlomann, Generation of phonons in high-power ferromagnetic resonance experiments, J.Appl.Phys. 31, 1647 (1960).
  * Eshbach [1963] J. Eshbach, Spin-wave propagation and the magnetoelectric interaction in yttrium iron garnet, J.Appl.Phys 34, 1298 (1963).
  * Schlomann [1964] R. Schlomann, E.and Joseph, Generation of spin waves in nonuniform magnetic fields.iii.magnetoelastic interaction, J.Appl.Phys 35, 2382 (1964).
  * Damon and van de Vaart [1965] R. W. Damon and H. van de Vaart, Dispersion of spin waves and magnetoelastic waves in yig, Proc. IEEE 53, 348 (1965).
  * Auld  _et al._ [1967] B. A. Auld, J. H. Collins, and H. R. Zapp, Adiabatic time domain conversion of hybrid magnetoelastic waves in yig, Appl. Phys. Lett. 10, 186 (1967).
  * Comstock and Kusnezov [1967] R. L. Comstock and N. Kusnezov, Magnetoelastic-elastic wave scattering, J. Appl. Phys. 38, 3740 (1967).
  * Rezende and Morgenthaler [1969] S. M. Rezende and F. R. Morgenthaler, Magnetoelastic waves in time-varying magnetic fields. i. theory, J. Appl. Phys. 40, 524 (1969).
  * Yukawa and Abe [1974] T. Yukawa and K. Abe, Fmr spectrum of magnetostatic waves in a normally magnetized yig disk, J. Appl. Phys. 45, 3146 (1974).
  * Camley [1979] R. E. Camley, Magnetoelastic waves in a ferromagnetic film on a nonmagnetic substrate, J. Appl. Phys. 50, 5272 (1979).
  * Komoriya and Thomas [1979] G. Komoriya and G. Thomas, Magnetoelastic-surface waves on yig substrate, J. Appl. Phys. 50, 6459 (1979).
  * Zaitsev  _et al._ [1998] B. D. Zaitsev, A. V. Ermolenko, and V. A. Fedorenko, Magnetoacoustic saw interaction in yig films, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 45, 356 (1998).
  * Filimonova  _et al._ [2004] Y. A. Filimonova, G. T. Kazakova, Y. V. Khivintseva, I. M. Kotelaynskiib, and A. V. Maryachinc, Nonlinear properties of magnetoelastic rayleigh waves in ferrite films, J. Magn. and Magn. Mater. 272, 1009 (2004).
  * Li  _et al._ [2008] M. Li, W. H. P. Pernice, C. Xiong, T. Baehr-Jones, M. Hochberg, and H. X. Tang, Harnessing optical forces in integrated photonic circuits, Nature 456, 480 (2008).
  * Weis  _et al._ [2010] S. Weis, R. Riviere, S. Deleglise, E. Gavartin, O. Arcizet, A. Schliesser, and T. J. Kippenberg, Optomechanically induced transparency, Science 330, 1520 (2010).
  * Teufel  _et al._ [2011] J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman, K. Cicak, A. J. Sirois, J. Whittaker, and S. R. Lehnert K.W., Sideband cooling of micromechanical motion to the quantum ground state, Nature 475, 359 (2011).
  * Chan  _et al._ [2011] J. Chan, T. P. Mayer Alegre, A. H. Safavi-Naeini, J. T. Hill, A. Krause, S. Groblacher, M. Aspelmeyer, and O. Painter, Laser cooling of a nanomechanical oscillator into its quantum ground state, Nature 478, 89 (2011).
  * Li and Zhu [2019] J. Li and S.-Y. Zhu, Entangling two magnon modes via magnetostrictive interaction, New J. Phys. 21, 085001 (2019).
  * Zhang  _et al._ [2020b] X. Zhang, G. Bauer, and T. Yu, Unidirectional pumping of phonons by magnetization dynamics, Phys. Rev. Lett. 125, 077203 (2020b).
  * Xu  _et al._ [2020] M. Xu, K. Yamamoto, J. Puebla, K. Baumgaertl, B. Rana, K. Miura, H. Takahashi, D. Grundler, S. Maekawa, and Y. Otani, Nonreciprocal surface acoustic wave propagation via magneto-rotation coupling, Sci Adv 6, eabb1724 (2020).
  * Vanderveken  _et al._ [2020] F. Vanderveken, F. Ciubotaru, and C. Adelmann, Magnetoelastic waves in thin films, arXiv preprint arXiv:2003.12099 (2020).
  * Ordonez-Romero  _et al._ [2007] C. L. Ordonez-Romero, O. V. Kolokoltsev, N. Qureshi, V. Grimalsky, and R. Ortega, Observation of indirect parallel pumping of magneto-elastic modes in layered yig/ggg structures, Solid State Commun. 141, 33 (2007).
  * Polzikova  _et al._ [2016] N. I. Polzikova, S. G. Alekseev, I. I. Pyataikin, I. M. Kotelyanskii, V. A. Luzanov, and A. P. Orlov, Acoustic spin pumping in magnetoelectric bulk acoustic wave resonator, AIP Adv. 6, 056306 (2016).
  * Chowdhury  _et al._ [2017] P. Chowdhury, A. Jander, and P. Dhagat, Nondegenerate parametric pumping of spin waves by acoustic waves, IEEE Magn. Lett. 8, 3108204 (2017).
  * An  _et al._ [2020] K. An, A. N. Litvinenko, R. Kohno, A. A. Fuad, V. V. Naletov, L. Vila, U. Ebels, G. de Loubens, H. Hurdequint, N. Beaulieu, _et al._ , Coherent long-range transfer of angular momentum between magnon kittel modes by phonons, Physical Review B 101, 060407 (2020).
  * Xu  _et al._ [2021] 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, Physical Review Applied 16, 024009 (2021).
  * Dutoit and Bellavance [1972] M. Dutoit and D. Bellavance, Ultrasonic attenuation measurements in gadolinium gallium garnet, in  _1972 Ultrasonics Symposium_ (IEEE, 1972) pp. 136–138.
  * Dutoit [1974] M. Dutoit, Microwave phonon attenuation in yttrium aluminum garnet and gadolinium gallium garnet, Journal of Applied Physics 45, 2836 (1974).
  * Chu  _et al._ [2017] Y. Chu, P. Kharel, W. H. Renninger, L. D. Burkhart, L. Frunzio, P. T. Rakich, and R. J. Schoelkopf, Quantum acoustics with superconducting qubits, Science 358, 199 (2017).
  * Multiphysics [2013] C. Multiphysics, Comsol multiphysics reference manual, COMSOL: Grenoble, France , 1084 (2013).
  * Gokhale  _et al._ [2021] V. J. Gokhale, B. P. Downey, D. S. Katzer, and D. J. Meyer, Phonon diffraction limited performance of fabry-pérot cavities in piezoelectric epi–hbars, in  _2021 IEEE 34th International Conference on Micro Electro Mechanical Systems (MEMS)_ (IEEE, 2021) pp. 206–209.
  * Chen  _et al._ [2019] H. Chen, N. F. Opondo, B. Jiang, E. R. MacQuarrie, R. S. Daveau, S. A. Bhave, and G. D. Fuchs, Engineering electron–phonon coupling of quantum defects to a semiconfocal acoustic resonator, Nano letters 19, 7021 (2019).
  * Schlitz  _et al._ [2022] R. Schlitz, L. Siegl, T. Sato, W. Yu, G. E. Bauer, H. Huebl, and S. T. Goennenwein, Magnetization dynamics affected by phonon pumping, arXiv preprint arXiv:2202.03331 (2022).
  * Feng and Winful [2001] S. Feng and H. G. Winful, Physical origin of the gouy phase shift, Optics letters 26, 485 (2001).
  * Renninger  _et al._ [2018] W. Renninger, P. Kharel, R. Behunin, and P. Rakich, Bulk crystalline optomechanics, Nature Physics 14, 601 (2018).
  * Fox and Li [1961] A. G. Fox and T. Li, Resonant modes in a maser interferometer, The Bell System Technical Journal 40, 453 (1961).
  * Poplavskiy  _et al._ [2018] M. V. Poplavskiy, A. B. Matsko, H. Yamamoto, and S. P. Vyatchanin, On fundamental diffraction limitation of finesse of a fabry–perot cavity, Journal of Optics 20, 075609 (2018).
  * Vinet and Hello [1993] J. Y. Vinet and P. Hello, Matrix simulation of optical cavities, Journal of Modern Optics 40, 1981 (1993).
  * Wang  _et al._ [2019b] Q. Wang, B. Heinz, R. Verba, M. Kewenig, P. Pirro, M. Schneider, T. Meyer, B. Lägel, C. Dubs, T. Brächer, _et al._ , Spin pinning and spin-wave dispersion in nanoscopic ferromagnetic waveguides, Physical review letters 122, 247202 (2019b).
  * Guslienko and Slavin [2005] K. Y. Guslienko and A. Slavin, Boundary conditions for magnetization in magnetic nanoelements, Physical Review B 72, 014463 (2005).
  * Kakazei  _et al._ [2012] G. Kakazei, G. Aranda, S. Bunyaev, V. Golub, E. Tartakovskaya, A. Chumak, A. Serga, B. Hillebrands, and K. Y. Guslienko, Probing dynamical magnetization pinning in circular dots as a function of the external magnetic field orientation, Physical Review B 86, 054419 (2012).
  * Matsushima and Shimobaba [2009] K. Matsushima and T. Shimobaba, Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields, Optics express 17, 19662 (2009).
  * Sokolov  _et al._ [2016] N. S. Sokolov, V. V. Fedorov, A. M. Korovin, S. M. Suturin, D. A. Baranov, S. V. Gastev, B. B. Krichevtsov, K. Maksimova, A. Grunin, V. Bursian, L. Lutsev, and M. Tabuchi, Thin yttrium iron garnet films grown by pulsed laser deposition: Crystal structure, static, and dynamic magnetic properties, Journal of Applied Physics 119, 023903 (2016).
  * Vitko  _et al._ [2015] V. V. Vitko, A. A. Nikitin, A. V. Kondrashov, A. A. Nikitin, A. B. Ustinov, P. Y. Belyavskiy, B. A. Kalinikos, and J. E. Butler, Microwave filter based on lamb modes for optoelectronic generator, in  _Journal of Physics: Conference Series_ , Vol. 661 (IOP Publishing, 2015) p. 012049.
  * Kharel  _et al._ [2018] P. Kharel, Y. Chu, M. Power, W. H. Renninger, R. J. Schoelkopf, and P. T. Rakich, Ultra-high-q phononic resonators on-chip at cryogenic temperatures, Apl Photonics 3, 066101 (2018).
  * Newberry and Thompson [1989] B. P. Newberry and R. B. Thompson, A paraxial theory for the propagation of ultrasonic beams in anisotropic solids, The Journal of The Acoustical Society of America 85, 2290 (1989).