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

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 and have diverse applications ranging from aerospace to medical imaging fields. Bulk acoustic wave (BAW) phonons are one of the most attractive candidates for quantum memories because of their long lifetime and high operating frequency. In this work, we establish a modeling approach that can be broadly used to design hybrid magnonic high-overtone bulk acoustic wave resonator (HBAR) structures for high-density, long-lasting quantum memories and efficient quantum transduction devices. We illustrate the approach by investigating a hybrid magnonic system, where BAW phonons are excited in a gadolinium iron garnet (GGG) thick film via coupling with magnons in a patterned yttrium iron garnet (YIG) thin film. We present theoretical and numerical analyses of the diffraction-limited BAW phonon lifetimes, modeshapes, and their coupling strengths to magnons in planar and confocal YIG/GGG HBAR structures. We utilize Fourier beam propagation and Hankel transform eigenvalue problem methods and discuss the effectiveness of the two methods to predict the HBAR phonons. We discuss strategies to improve the phonon lifetimes, since increased lifetimes have direct implications on the storage times of quantum states for quantum memory applications. We find that ultra-high, diffraction-limited, cooperativities and phonon lifetimes on the order of ~10^5 and ~10 milliseconds, respectively, could be achieved using a CHBAR structure with 10mum lateral YIG dimension. Additionally, the confocal HBAR structure will offer more than 100-fold improvement of integration density. A high integration density of on-chip memory or transduction centers is naturally desired for high-density memory or transduction devices.

The performance figure of merit for magnon-phonon transduction can be characterized using cooperativity, C = 4g 2  mb /κ m κ b , where g mb , κ m , and κ b are magnonphonon coupling strength, magnon dissipation rate, and phonon dissipation rate, respectively.The lifetimes of magnons and phonons are given by τ m = 1/κ m and τ b = 1/κ b , respectively.In YIG/GGG hybrid magnonic HBAR devices, the phonon lifetimes usually surpass the magnon lifetime thanks to the excellent mechanical properties of the garnets.Past studies reported the magnon and phonon lifetime in YIG/GGG HBAR systems to be τ m = 0.07µs and τ b ∼ 0.25µs, respectively [57,58].The system lifetime is primarily determined by the phonon lifetimes.It is desirable to further improve the phonon lifetimes and consequentially, the system lifetimes, for quantum memories and other information processing applications.However, a complete understanding about the loss mechanisms that limit the phonon lifetime in garnet devices is not yet available.The phonon lifetime could be limited by material and diffraction losses depending on the device geometry or the operating conditions.At room temperature, the phonon lifetime is primarily limited by acoustic attenuation due to phononphonon interactions [59,60].The acoustic attenuation effects, α ∝ 1/τ b , follows a T 4 temperature dependence.Due to the T 4 behavior, these effects are less significant at low temperatures, such as millikelvin regime where many qubit systems operate.Ideally, τ b can increase by multiple orders of magnitude at cryogenic temperatures.However, the diffraction losses play an important role in determining the phonon lifetime at low-temperature operating conditions, where the material losses are suppressed.In this study, we establish a computational approach for analyzing the diffraction losses in HBAR structures and use it to investigate the effects of diffraction losses on the HBAR phonons of hybrid YIG/GGG magnonic structures.We focus on the diffraction effects because they can be analyzed computationally, while experimental characterization is needed to investigate the acoustic attenuation effects.Note that the knowledge of material losses in YIG/GGG material systems is limited.There is a strong need for experimental characterization of acoustic attenuation in different magnetic materials across temperature and operating frequencies, to unlock their full potential.
The integration density of the hybrid magnonic HBAR structures is an important design aspect for developing high-density phononic quantum memories.The integration density is measured by the number of memory or transduction components that could be incorporated into a single chip, in addition to other on-chip circuitry.The YIG film represents the transduction component of the YIG/GGG HBAR structures of interest.Thus, the smaller the lateral dimension of the YIG film, the greater the number of transduction components in a single chip of given dimension and the integration density.We define the integration density as 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 for YIG/GGG devices in literature to be the reference, A 0 d = 0.8 × 0.9 mm 2 = 0.72 mm 2 [58].The integration density of the device can be increased by reducing the lateral dimension of the YIG films.However, the reduced aperture will lead to high-diffraction and affect phonon lifetime.In addition, a smaller aperture will reduce the overlap between the magnon and the phonon modes since they only overlap in the YIG film.A reduced magnon-phonon overlap could lead to a weaker magnon-phonon coupling strength.Consequently, the cooperativity will decrease as a result of the reduced lifetime and coupling strength.It is therefore imperative to develop a strategy to increase the integration density without inducing excessive diffraction losses.Note that we only focus on the lateral integration density while keeping the thickness fixed at 527.2µm for all structures investigated in this work.The same thickness allows us to keep the free spectral range of phonons fixed for all our analysis.
In this study, we investigate a set of planar and confocal YIG/GGG HBAR structures.We investigate the phonon lifetime and the magnon-phonon coupling strength in these structures, and identify strategies to improve their performance and integration density.We assume a magnon lifetime of 0.07µs [58] in all structures considered in this study.The phonon lifetime τ b or the quality factor (Q = ωτ b ) is of particular interest due to its direct relevance to the quantum memory storage time.Henceforth, we drop the subscript b from τ b , unless otherwise needed to distinguish it from magnon and photon lifetimes.Without losing generality, we only consider the fundamental magnon mode, i.e., the Kittel mode, for simplicity.The Kittel mode has a well-defined analytical function for circular discs.However, such analytical description is still lacking for the HBAR phonon modes.Several methods have been implemented to model HBAR phonons and their lifetimes.
For example, phonon lifetime in planar HBARs was calculated by decomposing the initial beam in a Bessel function basis and calculating the overlap of the reflected beam with the initial profile after each round trip [61].However, the predicted lifetime was smaller than the experimental observations because the effects of the lateral confinement induced by the transducer are not included.Finite element analysis (FEA) is another popular method to study HBAR phonons, using the COMSOL Multiphysics [62] software.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 [63].Another study used FEA to predict the modeshapes of a confocal HBAR for qubit coupling applications [64].They considered phonon modes at low-overtones or long wavelengths, hence it was possible to perform three-dimensional (3D) FEA simulations for this system.However, as we will discuss later, it becomes intractable to simulate bulk acoustic waves with overtone numbers n∼3000 using FEA since sampling required for these small wavelength modes is on the order of ∼ billion nodes.Moreover, all these studies discussed longitudinal HBAR phonons whereas in hybrid magnonic devices shear phonons exhibit much stronger coupling with magnon modes [57,58,65].The FEA analysis becomes particularly expensive for shear phonons since they are not axi-symmetric.One cannot leverage the axi-symmetric two-dimensional (2D) FEA modelling available in COMSOL Multiphysics for these systems.Due to these reasons, we explore other approaches to model shear phonon modes.We consider the Fourier beam propagation method (FBPM) [51,61] which follows a Fox-Li-like [66] 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 a detailed description of this method and discuss an adaptive algorithm that allows to overcome some of the challenges of the standard FBPM method.Additionally, we consider another method that uses Hankel transform (HT), which is a Fourier transfrom for axi-symmetric systems, and works with a Bessel function basis.Thus far, HT method has only been implemented for Fabry-Perot optical cavities [67,68].We adapt this approach to simulate YIG/GGG HBAR structures by leveraging their axi-symmetry and isotropic material properties.

HBAR Configurations
Figure 1 shows the two representative HBAR structures investigated in this work.The structure in Figure 1(a) consists of a thick GGG film joined 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.The thickness of the GGG film for all configurations considered in this study, is t GGG = 527µm, as shown in Fig. 1.All YIG films are circular films or discs with thickness, t YIG = 200 nm, and radius of cross-section, R YIG .The total HBAR device thickness is t HBAR = t GGG + t YIG = 527.2µm.We consider planar and confocal HBAR structures with several R YIG 's ranging from 10µm to 200µm.The width W of the structure is 1200µm unless otherwise mentioned.The radius of cross-section of the dome, R cross , and its radius of curvature, R curv , vary for different structures considered.We investigate 8 planar HBARs with varying R YIG , 18 CHBARs with varying R curv and fixed R cross , and 10 CHBARs with varying R cross and fixed R curv .
FIG. 1. Representative HBAR configurations: (a) Planar HBAR structure consisting of a GGG thick film and a YIG thin film.(b) Confocal HBAR structure with a top dome 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.Origin of the coordinate axes is at the center of the bottom YIG surface.

Magnon Modes
The YIG thin film hosts magnons, generated by a static, B 0 , and an oscillating RF magnetic field, B RF .B 0 is along the z axis, while B RF acts in the YIG plane.The magnetic fields, B 0 and B RF , generate forward volume magnetostatic modes in the YIG film, that precess in the x − y plane.The resulting dynamic magnetization is given by m YIG (x, y) = m 0 (x, y)(cos(ωt)i + sin(ωt)j).Here, m 0 (x, y) represents the magnon modeshape and ω is the frequency of precession.We practice the following convention throughout the article, bold-font symbols represent vector fields and corresponding normalfont symbols represent scalar values, respectively.We expect un-pinned spin waves at the top and bottom YIG film surfaces because of the following reason.The exchange interaction effects are expected to dominate at a thickness on the order of 200 nm and below.As a consequence, the pinning effect is expected to disappear below a critical width on the order of 200 nm [69].However, we expect full-pinning of the magnetic spin waves at the lateral boundaries of the YIG discs of radii ranging from 10 micrometers (µm) to 200µm.This is because the magnetic dipolar interaction effects dominate at length scales on the order of micrometers, over the exchange interactions [69].Taking these into account, we assume that m YIG is constant throughout the YIG film thickness.We describe the modeshape of the pinned spin waves, m 0 (x, y), using truncated Bessel functions [70,71]: where radial coordinate, R = x 2 + y 2 and J 0 is the 0 th order bessel function.The zeroes of ζ j with j = 0, 1, 2, ..., correspond to the fundamental and higher-order magnon modes respectively.In this article, we consider the fundamental magnon mode or the Kittel mode, whose amplitude is represented by the function, J 0 ( R RYIG ζ 0 ).The modeshape is constant along the z-direction.

Phonon Modes
The forward volume magnon modes, m 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 is determined by the precession direction of the Kittel mode.The shear displacements can be expressed as u YIG (x, y) = u 0 (x, y)(cos(ωt)i+sin(ωt)j), where u 0 (x, y) is the phonon modeshape and ω is the precession frequency.In the following, we refer to the shear displacements, u YIG (x, y), as 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 the reflected helical propagating waves interfere, they form rotating standing shear wave modes at specific frequencies.While the traveling chiral phonons are helical, the standing waves are not helical.The top GGG surface does not induce a π 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) HT eigenvalue problem approaches to analyze the shear phonon modes in the chosen HBAR configurations.
(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 [51,72].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 FBPM has been widely used to analyze beam propagation in the field of optics [72].The mathematical formulation of FBPM for acoustic waves is well established [51].Recently, it has been used to study phonons in planar [61] and CHBAR structures [51], adapting an iterative method similar to the Fox-Li approach [66].Here, we implement a reformulated iterative approach in which the propagation is calculated using projector and propagator operators.The reformulation allows us to achieve a seven-fold speed up of computation time.Our approach is particularly advantageous for isotropic systems such as YIG and GGG.
We assume an initial shear displacement field with modeshape The iterative procedure begins with this initial input field and the superscript represents iteration index (i = 1).The subscript in u (1) 0 (x, y) refers to the fact that the displacement is computed at z = 0. Note that the lateral phonon modeshape is the same as the magnon Kittel modeshape, as shown 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.The Fourier decomposition of the input beam is given by: Here, i is the iteration index and ũ(i) 0 is defined on a N × N (k x , k y ) grid which is the Fourier conjugate of the spatial N × N (x, y) grid.We calculate the projections of ũ(i) 0 (k x , k y ) along the different polarization directions, m = 1, 2, 3, using the projector, P m , and obtain the amplitudes, A (i) m , of the shear (m = 1, 2) and the longitudinal (m = 3) modes: with Here, dm is the polarization vector for the plane-waves propagating along (k x , k y , k z,m ), (j, k, l) refers to the Cartesian directions, (j, k, l) = (1, 2, 3) and Σ cyc and Π cyc are cyclic sum and product operators, respectively, that cycle through variables j, k, and l.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 HBAR , as shown in Fig. 1.The displacement field of the propagated beam at the upper GGG surface is given by with G m (k x , k y ) = dm (k x , k y )e ikz,m(kx,ky)tHBAR .(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 can be derived from their respective slowness surfaces [51] for each polarization and values of (k x , k y ).However, the functional relation can be simplified to Here, ω is the frequency of the initial wave and v m is the velocity of phonons with polarization m.We use the simplified relationship and the material properties of YIG and GGG, shown in Table .I, to compute k z,m for the propagating waves in our isotropic YIG/GGG HBAR structures.Since the properties of YIG and GGG are similar, 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.Note that we compute k z,m , A m (Eq.4), and u tHBAR (x, y) (Eq.5) only for the first iteration, for a given ω.This aspect results in the seven-fold computational speed-up mentioned earlier.We choose the longitudinal polarization d3 to be along the unit vector k and the shear polarizations, d1,2 , to be two mutually perpendicular unit vectors perpendicular to k.For the isotropic systems investigated in this article, d m 's are mutually perpendicular and A m reduces to 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 intro-duce absorbing boundary regions that completely attenuate any waves that enter these regions [61].In this study, we implement absorbing boundaries of thickness W Abs = 50 × λ, where λ is the wavelength of the input field.The blue shaded regions in Fig. 1 show the absorbing boundaries.To simulate the effect of absorbing boundaries, we multiply u (i) tHBAR (Eq.5) with a reflection operator R tHBAR , defined as W eff = W − 2W Abs is the effective width of the simulation window without the absorbing boundaries.We propagate the attenuated reflected wave further through a distance t HBAR , We implement the beam propagation by following the approach outlined in Eqs. 3 -5b, to complete a full round trip.After every round trip, we multiply the resulting complex displacement field by an additional phase R 0,m for each polarization m: 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 iteration.We repeat this process for N round trips.At the end of N round trips, we calculate the complex sum of the displacement fields, U 0 (x, y) = Σ N i=1 u (i) 0 (x, y) at z = 0. Using the interference sum, U 0 (x, y), we restart the iterative process with U 0 (x, y) as the initial beam of the next restart, u (1) 0 (x, y) = U 0 (x, y).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 due to destructive interference.We continue this process until the varation of the modeshape is within a chosen tolerance.We obtain a converged standing shear wave with displacement field, u 0 (x, y), as a final outcome.
In addition to the displacement profiles, we are interested in identifying the frequencies of the shear modes.These modes could couple with the Kittel magnon modes generated in the YIG thin-film (Eq.1), in a rotating coordinate system.The frequency overtones for plane-waves traveling along z-direction in the HBAR structures are expected to be ω m,n = 2π × nvm 2tHBAR , with n = 1, 2, 3, ...∞.Here, t HBAR is the thickness of the structure, velocity of wave is v m and m represents polarization.The overtones are separated by vm 2tHBAR , known as the free spectral range (FSR).However, unlike plane-waves, the diffracting waves traveling in the z direction do not have welldefined k z,m 's leading to Gouy phase effects [75].The diffraction results in the shift of resonance frequencies from the monochromatic plane-wave overtones.To identify the resonance frequencies, we select beam of frequencies from a chosen frequency window, propagate the beam in the structure, and calculate the intensities of the interference sums (I = A Re[U 0 (x, y)] 2 dA) at z = 0 after N round trips.The frequencies for which the intensities are maximum are the resonance frequencies of the shear waves in the HBAR structure.We focus on the shear modes with m = 2, however, we could have as well chosen m = 1 as the dominant shear polarization in the rotating coordinate system.We sweep through a frequency window of ω 0 ± 5 × FSR using 50 steps.Here, the frequency of interest, ω 0 = 2π × 9.825 GHz, corresponds to the frequency of the 2960 th overtone of a standing plane-wave in the HBAR structure.We choose this overtone to match with the structure investigated in a previous article [58].We consider a YIG/GGG HBAR structure with t YIG = 200 nm, R YIG = 200µm, t GGG = 527µm, and L x,GGG = L y,GGG = 1200µm, for this analysis.at frequencies near ω 0 .We obtain these displacement fields after 800 round trips.This simulation included one restart after 400 round trips.We perform multiple tests and check that this number is large enough such that the Gouy phase effects are not significant on the interference sums.We obtain the normalized displacement fields before the 1st and the 2nd restarts, U 1 0 (x, y) and U 2 0 (x, y), respectively, and compute the weighted deviation (WD) between them: |(U 1 0 (x, y) − U 2 0 (x, y))|/|Σ x,y U 1 0 (x, y)|.We monitor the WD between restarts to check for convergence.The color scale of Fig. 2 (d) shows that the 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 modeshapes 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 using a low-pass filter or a smoother phase change factor than the function, R 0,m (x, y) (Eq.7), used in this work.
To obtain the converged modes and identify the true resonant frequency near an intensity peak of interest, we further narrow the frequency sweeping range with finer sampling (∼ 1 kHz).We set the frequency to be ω HBAR 0 = ω 0 + δω, with δω = 2π × 2.157 kHz and restart the iterative process with Eq. 2 as the input beam, to identify the true modal frequency near ω 0 .If such narrowing is not done, the modeshape can change significantly after each restart due to the Gouy phase effects the beam incurs due to diffraction [75].These effects results in the detuning of overtones with respect to the plane-wave overtones.Figure 3 shows 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 ω 0 without performing the narrowed frequency sweeping discussed above.Both the real and imaginary parts of U (1) 0y and U (2) 0y show significant variations between the restarts.These results establish that the δω correction is necessary to obtain the converged modes.
However, the gradual narrowing of the frequency sweeping is a tedious process to identify the necessary fine sampling.We need to restart the iterative process multiple times especially 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 YIG as low as 10µm corresponding to a very low Fresnel number, N F = 1.06, indicating a high-diffraction regime.We investigate these structures to discuss how a reduced actuator lateral area affects  .
the phonon modes and their lifetimes.The reduced area promises to increase the integration density of the HBAR devices towards high-density memory applications.The high-diffraction regime of these structures introduces significant Gouy phase effects on the propagating waves and results in detuning of frequencies.It becomes challenging to identify a resonant frequency by merely narrowing the frequency sweeping range.We implement an adaptive FBPM algorithm to circumvent this challenge.Note that we implement absorbing boundaries in this study to avoid undesired boundary reflections.The phonon modeshapes can be well predicted using this implementation.However, their lifetimes can be dependent on the shape and size of the boundaries.We leave the extensive investigation of the effect of boundary conditions on phonon lifetimes for a future investigation.

Phonon Modeshape: Adaptive Fourier Beam Propagation
We formulate the FBPM approach described above as an eigenvalue problem: where RT is the round trip operator that includes all the operations described in the previous section (Eqs.3-7) and Λ is the eigenvalue corresponding to the eigenvector u 0 .It is possible to solve the eigenvalue problem for 2D problems (e.g., if HBARs are represented by planar 2D structures), however, it cannot be directly solved for 3D structures of our interest.We develop an iterative method to obtain the resonant mode that satisfies Eq. 8 with a real Λ.Λ is real for a standing wave mode.A complex Λ applies a global phase to the input modeshape after a round trip, and the phase prevents the waves to interfere constructively over multiple round trips.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. Start with a input beam, u , where 5. If R (i) ≤ tolerance (tol), end simulation and output final eigenvalue and modeshape of resonant modes, else continue to the next step.We choose tol = 10 −6 for all a-FBPM calculations of this study.
8. Repeat the process until step 5 is satisfied or the maximum number of N round trips is reached, at which point we set u 0 (x, y), and restart the iterative procedure.
Figure 4 illustrates the effectiveness of the a-FBPM algorithm in predicting the resonant modes of the HBAR structures.We show the results for two planar HBAR structures with R YIG = 200 µm (red) and 10 µm (blue), respectively.The residual, R, decays as a function of the number of steps, as we approach the resonant mode.The step-like features of the residual decay occur after every N R round trips 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.We calculate the complex interference sum after every N R round trips.We choose N R = 400 and 40 for the HBARs with R YIG = 200 µm and 10 µm, respectively.We choose a smaller number of round trips, N R, for the 10 µm case, since the HBAR with R YIG = 10 µm represents a highdiffraction structure and the modeshape decays rapidly for this case.As shown in Fig. 4, it takes a total of ∼650 and ∼1200 round trips (including restarts) to obtain converged resonant modes for the 10 µm and the 200 µm case, respectively.We continue the iterative process till R < 10 −8 to show stability of the resonant mode even after tolerance is reached.
(2) Phonon Modeshape Analysis: Hankel Transform Eigenvalue Problem Our a-FBPM approach avoids the frequency sweeping process and mostly overcomes the slow convergence issues of the standard FBPM approach.However, it is still challenging to obtain converged phonon modes for the confocal HBAR structures of interest, using the a-FBPM approach.Here, we discuss an alternate method that uses the Hankel transform (HT) approach, the axisymmetric equivalent of the Fourier transform method.The HT approach has been widely used in the field of optics, e.g., Fabry-Perot cavities [67].However, it has not been applied for acoustics problems, to the best of our knowledge.The HT method allows us to obtain converged phonon modes with significantly reduced computational cost.The approach leverages the axi-symmetry of the problem to reduce the 3D problem to a 2D one.For example, the 2D problem can be obtained by representing HBARs as planar 2D structures.The axi-symmetry conditions are applicable for the YIG/GGG HBAR structures due to their isotropic material properties and the circular cross-section of the YIG film.Note that the scalar field describing the dominant shear-component of u 0 is expected to obey axi-symmetry, however, the vector field of the shear acoustic phonon modes does not obey axi-symmetry.Here, we provide a brief description of the HT method.We encourage the interested reader to find a detailed description of the method elsewhere [67,68].
We cast the HBAR structures shown in Fig. 1 into axisymmetric representations about the z-axis.We transform the x − y coordinates to the radial coordinate r, with limit 0 ≤ r ≤ W/2.We write the Hankel eigenvalue problem as RT hk u(r) = Λu(r), where u(r) represents the axisymmetric scalar phonon modeshape, Λ is the corresponding eigenvalue and RT hk is the round trip operator.Since the dimensionality of the problem is reduced from 3D to 2D, the Hankel eigenvalue problem can be solved to obtain the phonon modes directly.We discretize r as where j = 1, 2, 3, .., N and ζ j (j = 1, 2, 3, .., N ) are the roots of the J 1 Bessel function.The round trip operator for the HT approach RT hk is given by where P, the N × N propagator matrix, is defined as: Here, H + is the HT whose matrix inverse is taken to be an approximation of the inverse HT, and G is the Green's function obtained from the Fourier transform of the Fresnel propagator in the paraxial approximation.R tHBAR , and R 0 are the N ×N diagonal reflection matrix operators at z = t HBAR and z = 0, respectively and the diagonal elements are defined as 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 axisymmetric constraints, it makes it possible to obtain the phonon modes and lifetimes of HBAR structures at a reduced cost.The expedited analysis allowed us to analyze a large class of HBAR structures, as we show in the Results and Discussion section.

Diffraction-Limited Phonon Lifetime
We estimate the diffraction-limited phonon lifetimes using three 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 with modeshape u is given by E ∝ |u| 2 dxdy.We only integrate over the dxdy element since the modeshape largely remains unaltered throughout the thickness.We assume that the elastic energy of the acoustic beam decays exponentially with the propagation time.For a cavity phonon mode with a total initial elastic energy E tot , we represent the elastic energy left in the cavity after one round trip as E in = E tot e −tRT/τ .Here, t RT (= 2t HBAR /v 2 ) is the time taken by shear waves to complete a round trip and τ is the lifetime of the phonon mode.On the other hand, the ratio of E in to E tot is given by E in /E tot = Λ 2 .Here, Λ is the real eigenvalue obtained using the a-FBPM simulations.Using this relation, we can obtain the lifetime of the phonon mode as Since these computations are performed on a finite x − y mesh grid, Λ, and consequently, τ can be sensitive to the mesh density.Hence, it is important to select an optimal grid to obtain the converged values of Λ.In Fig. 5, we show the variation of Λ with mesh density N x (= N y ) for both the low-diffraction (R YIG = 200µm) and the high-diffraction (R YIG = 10µm) cases.We find that the variation of Λ is much less than 1% when N x ≥ 1024.Consequentially, we choose N x = N y = 1024 to compute phonon modes and lifetimes for all cases considered here.
(B) Exponential curve fitting method: In this method, we calculate the phonon mode lifetime by evaluating the diffraction loss of the beam as it travels in the HBAR structure.We estimate the diffraction loss from the overlap of the propagated mode profile with the initial input profile, after each round trip.We estimate the overlap using the following function: Here, t is an integral multiple of the time taken for each round trip, t RT .We allow the initial beam to travel for multiple round trips 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 YIG = 200µm.The finite width of the YIG film (transducer) induces a localizing effect on the propagating beam.The localizing effect can be modeled by introducing a phase factor.The blue and the red lines shown in Fig. 6 correspond to the two cases when we obtained the overlap function with or without the consideration of the phase factor, respectively.The two different decays highlight the effect of the finite width of the transducer on the propagating beam.We find that I(t) decays at a much slower rate when the YIG-induced phase factor is considered, compared to the case without the phase factor.We fit the two decays with exponential functions and obtain the phonon lifetimes as τ wophase = 0.1 ms and τ wphase = 14.1 ms, respectively.The remarkable two-orders of magnitude difference between these two predicted lifetimes highlights the importance of including the transducer or actuator phase effects for such an analysis.Subsequently, we included the phase effects in all our analyses to obtain the phonon lifetimes.We find that the lifetime predicted with the phase effects, τ wphase = 14.1 ms, closely matches with τ of 15.5 ms, predicted using the eigenvalue method discussed earlier.A past study used the exponential curve fitting approach to estimate the lifetime of the longitudinal HBAR phonons in an AlN/Sapphire structure [61].However, they treated the phonons as superpositions of Bessel functions and did not account for the localization effects induced by the AlN transducer.(C) Clipping method: In this method, we obtain the phonon lifetime using a similar expression as used for the eigenvalue method, Eq. 13: τ = −2tHBAR v2ln(Ein/Etot) .Here, E tot is the total initial elastic energy for a cavity phonon mode and E in is the elastic energy left in the cavity after one round trip.We calculate E in , contained in a volume, V in , of a cylinder with radius of the YIG disc and the thickness of the HBAR structure, E in ∝ Vin |u| 2 dxdy.We consider a converged modeshape, obtained with a-FBPM or HT method.This method results in the lowest estimate of the phonon lifetimes since it assumes that all energy outside the cylinder of interest is lost after a round trip.Our approach is similar to the clipping method used for Fabry-Perot [67,68] and HBAR resonators, excited opto-mechanically using photon beams [76].In these optical systems, the clipping method is more applicable since the input photon beam spans the entire x − y 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.

Magnon-Phonon Coupling
The magnon-phonon coupling strength, g mb , is given by dz ) represent the complex conjugates of the x and y components of displacement u and their corresponding shear strains, respectively.The integration domain is the volume of the YIG film, V YIG , since the magnon modes reside in YIG.The term A ν Q η corresponds to the magnon and phonon normalizing constants.We normalize the magnon modes as Here, µ 0 M s = 0.175T is the saturation magnetization, γ = 2π × 28.5 GHz/T is the gyromagnetic ratio, k is the unit vector along the z axis, ν is the magnon mode index, and A ν is the magnon normalization constant, respectively.Similarly, we normalize the HBAR phonon modes as Here, ρ = 7080 kg/m 3 is the GGG material density, η is the phonon mode index, and ω η is the corresponding phonon frequency, respectively.V HBAR is the volume of the HBAR structure.Q η is the phonon normalization constant and is of the same dimensionality as A ν .
Note that if we assume zero diffraction, the acoustic energy is fully contained in the volume V in , the volume of the cylinder with radius of the YIG disc and the thickness of the HBAR structure.Such a scenario results in a complete overlap of lateral mode profiles of the magnon and the phonon modes.We obtain the zero diffraction limit of the magnon-phonon coupling strength to be g 0 mb /2π = 1.13 MHz.This value is independent of the width of the YIG film.Our result compares well with a previous experimental result of 1 MHz [57], and is slightly higher than another experimental result of 0.75 MHz [58].However, we consider the regime where diffraction causes the spread of acoustic energy outside V in , which can impact both τ and g mb of phonon modes in HBAR structures.We present the diffraction-limited HBAR phonon lifetimes and the magnon-phonon coupling strengths in the Results and Discussion section.

RESULTS AND DISCUSSION
We investigate hybrid magnonic HBAR structures that can support the development of high-density phononic quantum memories.We can include a greater number of transduction components in a single chip of given dimension by reducing the radius of the YIG film.However, the reduced aperture increases the diffraction of the acoustic waves propagating into the GGG region.Here, we discuss the diffraction-limited phonon lifetimes and the magnon-phonon coupling in the HBAR structures.

Diffraction-Limited Phonon Lifetime
Figure 7 (a) shows the variation of phonon lifetimes in planar HBAR structures as we decrease the YIG radius.We compute the lifetimes using the a-FBPM approach and following the eigenvalue method (red line), the exponential fitting (blue and purple lines), and the clipping methods (green line), respectively.In addition to the a-FBPM predictions, we show predictions from HT method in Fig. 7 using circles, with colors matching to their corresponding a-FBPM predictions.Note that we use |Λ| instead of Λ in Eq. 13 to calculate τ using the HT method.All methods predict that the lifetimes decrease with decreasing R YIG due to increased diffraction, as expected.The eigenvalue method results in the highest lifetime estimates, among the three methods considered.We obtain the highest lifetime to be τ = 15.5ms(Q = ωτ = 9.57 × 10 8 ) for HBAR structure with R YIG = 200µm.The lowest lifetime predicted by the eigenvalue method, is τ = 10.5µs(Q = 6.48 × 10 5 ), corresponding to HBAR structure with R YIG = 10µm, respectively.These results show that the lifetimes are reduced by more than three orders of magnitude in the high-diffraction regime.The lifetimes predicted by the eigenvalue method (red) closely match with those predicted by the exponential fitting method (blue), for the low-to-medium diffraction cases, when the phase effects are included.This result can be explained through the following argument.The input beam can be thought of 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 frequency 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 (exponential fitting method) is expected to have a similar character to that of the decay of the dominant eigenmode (eigenvalue method).As a result, we obtain similar lifetimes using the two methods.
However, in the high diffraction case (R YIG = 10µm), we observe 35-fold reduced lifetimes predicted by the exponential fit compared to that of the eigenvalue method.Due to high-diffraction, significant amount of the energy of the input acoustic beam spreads out laterally during the first few round trips, resulting in a rapid initial decay of I(t).In the exponential fit method, we obtain the lifetime from the exponential fit of I(t), starting at t = 0.The loss of acoustic overlap during the initial transient phase results in lower τ predictions.We expect that the two lifetime predictions will be closer if the exponential fit is obtained after a steady state is reached.Note that we obtain the exponential fitting predictions (blue) by including the localizing phase effects due to the YIG film.We also show the lifetimes predicted by this method, when we do not consider the phase effects due to the YIG film (purple).We find that the τ wophase 's are consistently lower than the τ wphase 's.The difference between the two predictions decreases monotonously from 135 to 1.35 fold, as we decrease R YIG from = 200µm to = 10µm, respectively.This is expected since the localizing effect of the YIG film diminishes as we approach R YIG → 0. Past studies often ignored the effect of aperture width on the HBAR phonons, however, our results (Fig. 6 and Fig. 7) establish that such effects needs to be considered to make reliable predictions.
Figure 7 (a) also shows 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 obtain the lowest lifetime in the highdiffraction regime to be τ = 0.311µs (Q = 1.92 × 10 4 ).Interestingly, this prediction is still marginally greater than the experimentally observed values of 0.25µs and 1.45 × 10 4 , respectively [58], for a device operating in the low-diffraction regime.The reason for this discrepancy is that the experimental values correspond to HBAR devices operating at room temperatures.In this limit, the performance of HBAR structures is limited by the material attenuation effects [59,60].In our study, we ignore the material attenuation effects and only discuss the diffraction-limited performance.Our study will correspond to devices operating in the cryogenic or potentially milli Kelvin regimes.For example, our results will have high applicability for the HBAR devices that could effectively couple with superconducting qubit systems, etc., which typically operate in the milli Kelvin regime.
We discuss 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 HBAR ) and a function of the phonon velocity and the eigenvalue, Λ, as shown in Eq. 13.Following Eq. 13, an FIG. 7. Diffraction-limited properties of phonon modes of planar HBAR structures: (a) Phonon lifetimes, τ , 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).Reduction of g mb as we decrease RYIG, is connected to the spread of modeshape in the high-diffraction regime.Increasing amount of acoustic energy spreads into the GGG region, causing reduced overlap between phonon and magnon modes in the YIG region.The ratio of acoustic energy spreading out to the total energy, Eout/Etot (black), increases in the high-diffraction regime.(c) Magnon-phonon cooperativity C. RYIG 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 HT predictions.
error propagation relation can be written as The uncertainty dτ τ increases monotonously with τ when other factors are fixed.It can be deduced from Eq. 18 that to predict a lifetime within x% of variability, Λ has to be predicted to be within tHBARx v2τ % of variability.Here, the order of magnitude of the prefactor v2 tHBAR ∼ 10 7 .The prefactor remains fixed since we keep the material and the HBAR thickness unaltered.This implies that we have to be increasingly more stringent with our Λ convergence criteria for high lifetime calculations.Figure 8 shows the conditions for desired accuracy (relative tolerance) needed for Λ predictions for high-τ (high-Q) systems.When τ = 0.1µs and the desired accuracy is set to 10%, Λ 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 τ = 1ms, and the desired accuracy is set to 10% of τ value, Λ must be predicted within ∼10 −3 % of accuracy.For our R YIG = 200µm HBAR structure, we continue the simulation till we obtain dτ τ = 7.9 × 10 −4 with dΛ Λ = 4.7 × 10 −7 , between the last two restarts.On the other hand for the HBAR with R YIG = 10µm, we continue till dτ τ = 3.9 × 10 −7 with dΛ Λ = 3.8 × 10 −8 .However, achieving such high accuracy requires high computational expense associated with simulating a large number of iterations.We implement and use an alternative HT approach to expedite the numerical analysis.

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 estimate the effect of diffraction by computing the ratio between the coupling strength in the planar HBAR structures, g mb , and that in the zero-diffraction limit, g 0 mb .In Fig. 7(b)(red), we show the variation of the ratio, g mb /g 0 mb , with decreasing YIG radius.As expected, The ratio is equal to 1 for low-diffraction cases with R YIG ≳ 100µm.As we reduce R YIG to 40µm, g mb only changes by 5% from the g 0 mb limit.Note that even though this case typically corresponds to a strong-diffraction regime, the effect of diffraction on g mb is minimal.This is because the YIG disc helps to localize the beam and thus, preserve the magnon-phonon overlap and the coupling strength.However, as we decrease R YIG ≤ 40µm, we observe a sharp decline of g mb .The percentage of the ratio, g mb /g 0 mb , drops to 77.4% and 55.7% for HBAR with R YIG = 20µm and 10µm, respectively.The decrease of g mb can be explained in the following way.In the zero diffraction limit, the acoustic beam is completely localized in the volume V in , of a cylinder with radius equal to R YIG and the thickness of the HBAR structure.This results in a complete overlap of lateral mode profiles of magnon and phonon modes.Also, the acoustic energy of a propagating beam is completely contained in the volume V in , E tot = E in .However, as we decrease R YIG , the modeshape spreads out beyond V in , due to diffraction.As a result, there is reduced overlap between the phonon and the magnon mode profiles in the YIG region.Correspondingly, some acoustic energy also leaks out of V in and spreads into the GGG region.We refer to the energy of the beam lying outside V in as E out . Figure 7(b) (black) shows that the ratio E out /E tot increases in the high-diffraction regime, with decreasing R YIG .
We combine the magnon-phonon coupling strength (Fig. 7(b)) and the phonon lifetimes (Fig. 7(a)) to determine the performance figure of merit of the HBAR structures, defined by the magnon-phonon cooperativity, We show the variation of C with R YIG in Fig. 7(c).We assume the magnon lifetime to be τ m = 0.07µs [57,58].We compute the cooperativity values using τ b predicted from eigenvalue method (red line), exponential fitting (blue and purple lines), and clipping methods (green line), respectively.We obtain a monotonically decreasing diffraction-limited C ranging from 21.9×10 4 to 46.4, using τ from the eigenvalue method, as we decrease R YIG from 200µm to 10µm.On the other hand, C is in the range between 1.2 × 10 4 and 1.4, predicted using τ from the clipping method.As can be noted from Fig. 7, τ , g mb , and C predicted by the HT method are in an excellent agreement with a-FBPM predictions, for all R YIG cases considered.The close match validates our implementation of the HT method.As we have discussed earlier, a-FBPM has broader applicability compared to HT method, and can be applied for systems with anisotropic material properties and transducer shapes.However, we find that it is an expensive and uncertain process to converge the residual since 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, HT 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 HT method in this article.
Note that our calculations predict high cooperativity for even the HBAR with R YIG = 10µm.This implies that if material acoustic attenuation effects are eliminated (e.g. by operating in the milli Kelvin regime), one can achieve high cooperativity for planar HBAR structures.However, the phonon lifetime of the HBAR structure with R YIG = 10µm is limited to 10.5µs (0.31µs), as predicted by the eigenvalue (clipping) methods and the maximum magnon-phonon coupling strength is only 55.7% of g 0 mb .We explore design approaches to further improve these performance parameters.
Enhancing integration density and Performance in Diffraction-Limited Regime FIG. 9.
Variation of magnon-phonon coupling, g mb /g 0 mb , with varying waist of fundamental phonon mode: g mb /g 0 mb is maximum when w0/RYIG ∼ 0.65.
Here, we illustrate that the use of focusing dome-like surface structures could significantly improve the performance of HBAR structures.Past studies demonstrated that the confocal HBAR structures could achieve Qfactors on the order of 10 7 , while resulting in 10 3 -fold reduction in device volumes [76].However, such structures have not been explored for hybrid magnomechanical systems.We show that both phonon lifetimes and magnon-phonon coupling can be improved by employing confocal geometries, leading to improved performance of hybrid magnomechanical systems.We first identify the optimal domeshape that will result in the maximum overlap between the CHBAR phonon and the magnon modeshapes, and consequently, the highest g mb .We obtain theoretical estimates of the overlap and identify the dome shape where the overlap is maximum.We assume that the phonon modeshapes in CHBAR structures can be described by Laguerre-Gaussian (LG) functions.The fundamental mode, LG 00 , can be assumed to be of the FIG.10.Performance of CHBAR structures with varied radius of curvature, Rcurv, of the dome-shape: (a) Phonon lifetimes, τ , calculated from the eigenvalue (red) and clipping methods (green), following the HT method, (b) Ratio between magnon-phonon coupling in CHBAR structures and that in the zero-diffraction limit, g mb /g 0 mb and (c) Magnon-phonon cooperativity, C. Rcurv is varied keeping RYIG = 10µm and Rcross = 60µm fixed.Solid lines represent corresponding values for a planar HBAR structure with RYIG = 10µm, while the circles of the same color correspond to the CHBAR form: u y ∝ e −R 2 /w 2 0 , due to the axi-symmetry of our system.Here, R is the radial coordinate x 2 + y 2 .w 0 is the waist of the Gaussian beam at z = 0 (see Fig. 1), that can be varied by changing t HBAR , ω, or the dome radius of curvature, R curv .We normalize the modeshapes according to Eqs. 16, 17 and calculate the coupling strengths using Eq. 15. Figure 9 shows the computed g mb /g 0 mb ratio as a function of the beam waist, w 0 .We find that the peak of g mb /g 0 mb occurs at w 0 /R YIG ∼ 0.65.g mb decreases sharply if w 0 ≲ 0.65R YIG , however, decreases slowly if w 0 is increased beyond the optimal value.Interestingly, we find that w 0 /R YIG is a constant and does not depend on R YIG , for all the CHBAR structures considered in this article.We use the w 0 value to obtain an initial estimate of the radius of curvature, R curv , of the dome [77]: with We obtain the analytical estimate to be R curv /t HBAR ∼1.5 that corresponds to w 0 /R YIG ∼0.65 and the peak of g mb /g 0 mb as shown in Fig. 9.Note that R curv can be similarly estimated for various device thicknesses (t HBAR ), wavelengths (λ) and waists by appropriately accounting for these variables shown in Eq. 19.This R curv estimate results in maximum g mb /g 0 mb , however, this analysis does not consider the effects of the dome shape on phonon lifetimes and cooperativities.Also, we ignore the effects of the lateral extent of the HBAR structure, the localizing effects of the YIG film, and assume phonon modeshapes to be perfectly Gaussian.
To obtain a R curv estimate that improves the overall performance of these CHBAR structures, we carry out a systematic numerical analysis using the HT method.
We modify the reflection operator at the upper GGG surface to include the effects of the lateral extents of the dome and the device, the attenuation window, and the localization effects of the YIG film.We include the phase induced by the dome surface in Eq. 6 and the radial form of Eq. 11c and write the modified reflection opertor as We vary R curv /t HBAR across the analytical estimate, ∼1.5, and identify R curv that results in the optimal performance of the CHBAR structures.We show the effect of the focusing dome on the phonon lifetime, magnonphonon coupling and cooperativity of the CHBAR structure in Fig. 10(a), (b) and (c), respectively.We show the predictions from the eigenvalue method and the clipping method using red and green circles, respectively.We also show the corresponding performance parameters of the planar HBAR structure with red and green horizontal lines, respectively.For the clipping method, we consider the volume, V in , defined by the dome's lateral cross-sectional area and the CHBAR thickness.We fix R YIG = 10µm and the dome x − y cross-section radius, R cross = 60µm, for this analysis.As can be seen from Fig. 10(a), the phonon lifetimes are maximum at R curv /t HBAR = 1.09.This value is less than the analytical prediction for maximum g mb , R curv /t HBAR ∼ 1.5.However, the peak value, τ ∼218ms (red circles), is ∼500 times larger than the value at R curv /t HBAR ∼1.5, justifying the necessity of performing the systematic analysis.The τ peak value for the CHBAR structure is ∼1.7 × 10 4 (red circles) (∼7.3 × 10 4 , green circles) times greater than the corresponding planar HBAR τ prediction, shown with red (green) horizontal lines, respec- tively.The peak τ predicted by the clipping method is ∼22.7ms,lower than those predicted by the eigenvalue method, as expected.We note that τ decrease sharply for R curv /t HBAR <1.This corresponds to the situation when the center of curvature of the dome is inside the HBAR structure, resulting in a negative 'g−parameter' [76], where g = 1 − tHBAR Rcurv .It is shown that one cannot obtain real and finite Gaussian beam solutions for structures with a negative g−parameter [76].We find that this is the case for the modes corresponding to R curv /t HBAR <1 in our case, that is, they are lossy modes which do not have well-defined Gaussian characters.
Figure 10(b) shows that the coupling, g mb /g 0 mb , is approximately 0.99 at R curv /t HBAR ∼1.4, a value close to the analytical estimate of ∼1.5.These results show that reducing the phonon mode volume by means of focusing does not necessarily improve g mb beyond the maximum achievable value, g 0 mb , which represents the case when there is complete lateral overlap of phonon and magnon modes.We use τ and g mb to calculate the figure of merit, C, and show the results in Fig. 10(c).The peak C value is predicted to be 2.5 × 10 6 (2.6 × 10 5 ), using the eigenvalue (clipping) method.The C values reach peak at R curv /t HBAR = 1.09,where τ is maximum, as we have noted from Fig. 10(a).The τ and C values of the CHBAR structure are several orders of magnitude improved compared to its planar counterpart and the coupling reaches up to 90% of g 0 mb .This result highlights that τ is the dominating factor that controls the performance of CHBAR structures, since g mb is always limited to the maximum value, g 0 mb .Note that we only vary R curv for this analysis and keep the radius of the dome cross-section, R cross , fixed at 60µm.Due to the fixed cross-sectional area of the dome, the lateral integration density of these devices, Here, A d is the cross-sectional area of the confocal dome surface and A 0 d = 0.72 mm 2 [58].To further improve the integration density, we analyze CHBAR structures with reduced dome cross-sectional area.We keep R curv /t HBAR = 1.09 and R YIG = 10µm fixed for this analysis.Figure 11 (a), (b) and (c) show the variation of τ , g mb , and C with R cross /R YIG , respectively.As we decrease R cross /R YIG from 6 to 1, τ is reduced by more than 4 orders of magnitude from ∼218 ms (22.7 ms) to 4.83µs (0.134µs), as predicted by the eigenvalue (clipping) method.g mb shows a sharp decline for R cross /R YIG <3 reaching the lowest value of g mb /g 0 mb ∼0.32 at R cross /R YIG = 1, however, remains mostly unchanged for R cross /R YIG ≥ 3. The cooperativity follows a similar trend to that of τ , 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×10 6 (2.6×10 5 ) to 7.2 (0.19) as we decrease R cross /R YIG .Additionally, τ , g mb , and C of these CHBAR structures is reduced below those of the planar HBAR values, shown with solid lines, if R cross /R YIG <2.Considering all the different aspects, we argue that the optimal geometry for the CHBAR structures is represented by the condition R cross /R YIG = 4.5.Our analysis using the HT eigenvalue (clipping) method predicts that the diffraction-limited lifetime of τ =144.7 ms (11.1 ms) and a integration density of D i = 113 could be achieved at a cooperativity C = 1.64×10 6 (1.25×10 5 ) for the CHBAR structure with R curv /t HBAR = 1.09,R YIG = 10µm and R cross /R YIG = 4.5.The results presented in Figs. 10 and 11 provide a proof-of-concept that including a focusing dome improves the performance of the hybrid magnonic HBAR structures for quantum memory and transduction applications.

CONCLUSION AND OUTLOOK
In this article, we establish a modeling approach that can be broadly used to design hybrid magnonic HBAR structures for high-density, long-lasting quantum memories and efficient quantum transduction devices.We illustrate the approach by discussing the magnon-phonon transduction properties of hybrid magnonic YIG/GGG HBAR structures.We present analytical and numerical analyses of the bulk acoustic wave phonon mode lifetimes, and the magnon-phonon coupling strengths in planar and confocal YIG/GGG HBAR structures.We discuss strategies to improve the phonon mode lifetimes of these structures, since increased lifetimes have direct implications on the storage times of quantum states for quantum memory applications.Additionally, high integration density of on-chip memory or transduction centers is naturally desired for high-density memory or transduction devices.In our structures, the transduction centers are represented by the YIG films and thus, integration density can be increased by reducing the lateral dimension of these films.However, the reduced aperture results in high-diffraction that affects the performance of these devices.We analyze the diffraction-affected shear wave phonon modes in the HBAR structures using two different methods: (1) Fourier beam propagation (FBPM) and (2) Hankel transform (HT) eigenvalue problem method.The FBPM method 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.We find that the FBPM method often requires us to analyze large and unpredictable number of round trips to obtain converged phonon modes.We implement an adaptive FBPM (a-FBPM) approach that mostly overcomes the slow convergence issues.However, we find that a-FBPM still suffers from convergence challenges when applied to confocal HBAR structures.To circumvent the challenges, we implement the HT method that allows us to obtain converged phonon modes with significantly reduced computational cost.The HT method leverages the axisymmetry of the problem and the isotropic material properties of the YIG/GGG HBAR structures to reduce the 3D problem to a 2D one and expedite the analysis.The HT method has been mostly used in the field of optics, e.g., Fabry-Perot cavities, however, it has not been applied for acoustics analysis, to the best of our knowledge.
Our study provides key insights into the diffractionlimited performance of the YIG/GGG HBAR structures.We predict the diffraction-limited τ to be on the order of milliseconds, for a planar HBAR structure with lateral YIG dimension, R YIG = 200µm.A recent study reported that the performance of a YIG/GGG HBAR structure at room temperature is limited by the phonon lifetime at 0.25µs [58].The previously studied structure had larger YIG lateral area than our structures and thus, the diffraction effects were less dominant.The phonon lifetime in HBAR structures could be limited by both material and diffraction losses depending on the device geometry and/or the operating conditions.At room temperature, the phonon lifetime is primarily limited by acoustic attenuation effects.However, these effects are less significant at low temperatures while the diffraction losses play an important role.Therefore, our results will have high applicability for devices operating in the cryogenic or potentially milliKelvin (mK) regimes.For example, our approach and analyses can be applied to design HBAR devices that could effectively couple with superconducting qubit systems.We acknowledge that the analysis of the material losses is necessary to obtain a complete understanding of the performance of the YIG/GGG systems.This may include an understanding of the different attenuation effects (e.g., Akhiezer and Landau-Rumer) at various temperatures and frequencies of interest, and strategies to overcome them.
Assuming that the material-limited lifetimes could be increased to ∼0.1 ms at mK temperatures, we find that the planar HBAR structures are not affected by diffraction effects even at R YIG = 50µm,.This structure already offers a significant 50-fold improved integration density over the reference structure [58].The integration density can be further improved by scaling down the YIG film lateral area.However, the reduced aperture will affect the phonon lifetime and the magnon-phonon coupling strength, resulting in a decrease of the cooperativity, the performance figure of merit for magnonphonon transduction.It is therefore imperative to develop a strategy that increases the integration density of HBAR structures without affecting the performance.Additionally, the performance of planar HBAR structures was found to be highly sensitive to the parallel nature of the surfaces.To address both these aspects, we investigate confocal YIG/GGG HBAR structures with top focusing domes and a planar bottom surfaces.Use of focusing dome structures has been proposed to significantly improve the phonon lifetime and integration density of these systems.Furthermore, the dome structures eliminate the necessity of maintaining perfectly parallel surfaces of planar HBAR structures.We first theoretically estimate the shape of the dome structure for which the magnon-phonon coupling strength, g mb , is at its maximum value.We then perform a rigorous numerical analysis and obtain a refined set of shape parameters of the dome for which both τ and C are optimal, and g mb is close to its peak value.Overall, we find that ultra-high, diffraction-limited, cooperativities and phonon lifetimes on the order of ∼10 5 and ∼10 ms, respectively, could be achieved using a CHBAR structure with R YIG = 10µm.In addition to enhanced τ and C, the confocal HBAR structure will offer more than 100-fold improvement of integration density.
In this work, we discuss the diffraction-limited performance induced by the lateral area of the YIG film and the dome structure.It will be interesting to apply the insights presented in the article and explore the applicability of YIG/GGG HBAR structures for quantum transduction applications, as future work.For example, these HBAR structures could be used for coupling with other quantum information carriers such as superconducting qubits, which operate in the microwave frequency regime and at mK temperatures.However, to optimize the performance of such devices, a further comprehensive analysis of different geometric and physical parameters will be necessary.For example, we keep the thickness fixed at 527.2µm for all YIG/GGG HBAR structures investigated in this work.The same thickness allows us to keep the free spectral range of phonons fixed for all our analysis.The different geometric parameters, such as the overall thickness, the thickness of the YIG and the GGG regions, geometrical misalignment, and imperfections, or, the physical parameters, such as the photon and magnon lifetimes and their coupling strengths, could also play important role in the overall device performance.We assume that these parameters remain invariant in our analysis.A comprehensive optimization analysis can be performed using the current state-of-the-art machine learning techniques, which is a promising research direction for the future.Additionally, we only discuss the coupling between the fundamental magnon mode and the fundamental phonon modes of the planar and the confocal HBAR structures.It will be interesting to consider the effect of higher-order mode couplings in the device performance.The higher-order mode couplings can be particularly relevant here since these structures host long-lasting phonon modes with broadband nature.Other interesting directions to explore are the anharmonic phonon interactions at the surface and interfaces, and the coupling of bulk acoustic waves and the surface acoustic waves.Our study is concerned with coherent states, and it will be important to explore how the insights provided here translate to systems involving non-classical states, e.g.squeezed states, cat states, and Fock states.

FIG. 2 .
FIG. 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 U0(x, y), U0y(x, y).(c).Top view of U0y(x, y) showing localization effect caused by the YIG film.(d) Weighted deviation (WD) between mode profiles before first and second restarts.

Figure 2 (
Figure 2(a) shows the intensities of different propagated beams with frequencies within the frequency window, a frequency window, ω 0 ± 5 × FSR.The central intensity peak corresponds to ω 0 = 2π × 2960v2 2tHBAR , the exact value of the 2960 th plane-wave overtone.The peaks form at resonance frequencies separated by FSR = 3.32 MHz around the frequency, ω 0 /2π ∼ 9.825 GHz, indicating

FIG. 3 .
FIG. 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 U (1) 0y and second restarts U at ω0.The components deviate significantly between restarts.
, y), for the i th round trip.The initial input beam (i = 1) has wavelength, λ (1) = 2t n and frequency ω (1) = 2π × nv2 2t .Here, n is the overtone number, t = t HBAR and v 2 is the velocity of shear phonons.The wavelength and frequency are estimated based on the n th overtone for an open-ended column.The inital modeshape is chosen to be a truncated Bessel function as shown in Eq. 2. 2. Calculate u (i+1) 0 (x, y) = RTu (i) 0 (x, y).

FIG. 4 .
FIG. 4. Obtaining resonant modes using the adaptive FBPM algorithm: Residual (R) of modeshapes decays as a function of number of round trips.Results are shown for two planar HBAR structures with RYIG = 200 µm (red, left and bottom axes) and 10 µm (blue, right and top axes).Circles with arrows point to corresponding axes for each plot.

FIG. 5 .
FIG. 5. Eigenvalues of two planar HBAR structures with R YIG = 200µm (red) and 10µm (blue): Variation of Λ for meshes with density ranging from 128 × 128 to 2400 × 2400.The width W is kept fixed at 1200µm.The circles with arrows point to the axes corresponding to each plot.

FIG. 6 .
FIG. 6. Decay of the overlap function, I(t), for propagating beams in a HBAR structure with R YIG = 200µ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).

= 7 ×
10 5 J/m 3 , is the magnetoelastic constant of YIG, (m x , m y ) are the x and y components of the magnetization, m = m YIG , (u * x , u * y ) and (

FIG. 11 .
FIG. 11.Performance of CHBAR structures with varied radius of cross-section, Rcross, of the dome-shape: (a) Phonon lifetimes, τ , calculated from the eigenvalue (red) and clipping methods (green), following the HT method, (b) Ratio between magnon-phonon coupling in HBAR structures and that in the zero-diffraction limit, g mb /g 0 mb and (c) Magnon-phonon cooperativity, C. Rcross is varied keeping RYIG = 10µm and Rcurv/tHBAR = 1.09 fixed.Solid lines represent corresponding values for a planar HBAR structure with RYIG = 10µm, while the circles of the same color correspond to the CHBAR values.

TABLE I .
Material properties of YIG and GGG.