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

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

$$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)

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:

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

(3)

Here, $i$ is the iteration index and $\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. We calculate the projections of $\tilde{\textbf{u}}^{(i)}_{0}(k_{x},k_{y})$ along the different polarization directions, $m=1,2,3$, using the projector, $\textbf{P}_{m}$, and 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}(j-m)|}\hat{d}_{j}(\hat{d}_{m}\cdot\hat{d}_{j}-\hat{d}_{k}\cdot\hat{d}_{l})}{1+2\Pi_{cyc}(\hat{d}_{j}\cdot\hat{d}_{k})-\Sigma_{cyc}(\hat{d}_{j}\cdot\hat{d}_{k})^{2}}.$		(4b)

Here, $\hat{d}_{m}$ 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 $\Sigma_{cyc}$ and $\Pi_{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_{\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}$ 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 $k_{z,m}=\sqrt{\frac{\omega^{2}}{v^{2}_{m}}-k^{2}_{x}-k^{2}_{y}}$, for isotropic dispersions. Here, $\omega$ 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. 1, 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 $\textbf{u}_{t_{\text{HBAR}}}(x,y)$ (Eq. 5) only for the first iteration, for a given $\omega$. This aspect results in the seven-fold computational speed-up mentioned earlier.

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. For the isotropic systems investigated in this article, $d_{m}$’s are mutually perpendicular and $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 [61]. In this study, we implement absorbing boundaries of thickness $W_{\text{Abs}}=50\times\lambda$, where $\lambda$ 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 $\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)

$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}}$, 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$:

$$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 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, $\textbf{U}_{0}(x,y)=\Sigma^{N}_{i=1}\textbf{u}^{(i)}_{0}(x,y)$ at $z=0$. Using the interference sum, $\textbf{U}_{0}(x,y)$, we restart the iterative process with $\textbf{U}_{0}(x,y)$ as the initial beam of the next restart, $\textbf{u}^{(1)}_{0}(x,y)=\textbf{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, $\textbf{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 $\omega_{m,n}=2\pi\times\frac{nv_{m}}{2t_{\text{HBAR}}}$, with $n=1,2,3,...\infty$. Here, $t_{\text{HBAR}}$ is the thickness of the structure, 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, unlike plane-waves, the diffracting waves traveling in the $z$ direction do not have well-defined $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=\int_{A}\text{Re}[\textbf{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 $\omega_{0}\pm 5\times\text{FSR}$ using 50 steps. Here, the frequency of interest, $\omega_{0}=2\pi\times 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_{\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(a) shows the intensities of different propagated beams with frequencies within the frequency window, a frequency window, $\omega_{0}\pm 5\times\text{FSR}$. 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}}$ plane-wave overtone. The peaks form at resonance frequencies separated by FSR $=3.32$ MHz around the frequency, $\omega_{0}/2\pi\sim 9.825$ GHz, indicating the presence of multiple shear wave overtones. In Figs. 2 (b) and (c), we show different views of the $y$-component of the normalized resonant modeshape, Re[U${}_{0y}(x,y)$], at frequencies near $\omega_{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, $\textbf{U}^{1}_{0}(x,y)$ and $\textbf{U}^{2}_{0}(x,y)$, respectively, and compute the weighted deviation (WD) between them: $|(\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. 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 ($\sim 1$ kHz). 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, to identify the true modal frequency near $\omega_{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 $\omega_{0}$ without performing the narrowed frequency sweeping discussed above. Both the real and imaginary parts of $\text{U}^{(1)}_{0y}$ and $\text{U}^{(2)}_{0y}$ show significant variations between the restarts. These results establish that the $\delta\omega$ 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_{\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 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.

We formulate the FBPM approach described above as an eigenvalue problem:

$$\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 operations described in the previous section (Eqs. 3- 7) and $\Lambda$ is the eigenvalue corresponding to the eigenvector $\textbf{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 $\Lambda$. $\Lambda$ is real for a standing wave mode. A complex $\Lambda$ 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. 1.

   Start with a input beam, $\textbf{u}^{(i)}_{\text{0}}(x,y)$, for the $i^{th}$ round trip. The initial input beam $(i=1)$ has 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 the $n^{\textbf{{th}}}$ overtone for an open-ended column. The inital 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 from, $\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⁽ⁱ⁾ $\leq\text{tolerance (tol)}$, end simulation and output final eigenvalue and modeshape of resonant modes, else continue to the next step. We choose $\text{tol}=10^{-6}$ for all a-FBPM calculations of this study.

6. 6.

   If R${}^{(i)}>\text{tol}$, estimate the global phase factor introduced 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)}=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 $N$ round trips 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 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_{\text{YIG}}=$ 200 $\mu$m (red) and 10 $\mu$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 $NR$ 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 $NR$ round trips. 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 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. We continue the iterative process till R $<10^{-8}$ to show stability of the resonant mode even after tolerance is reached.

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 axi-symmetric 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 $\textbf{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 axi-symmetric representations about the $z$-axis. We transform the $x-y$ coordinates to the radial coordinate $r$, with limit $0\leq r\leq W/2$. We write the Hankel eigenvalue problem as $\textbf{RT}_{\text{hk}}u(r)=\Lambda u(r)$, where $u(r)$ represents the axisymmetric scalar phonon modeshape, $\Lambda$ is the corresponding eigenvalue and $\textbf{RT}_{\text{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

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

(9)

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

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

(10)

where P, the $N\times N$ propagator matrix, is 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 HT whose matrix inverse is taken to be an approximation of the inverse HT, 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 the $N\times N$ diagonal reflection matrix operators at $z=t_{\text{HBAR}}$ and $z=0$, respectively and the diagonal elements 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)

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 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.

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\propto\int|\textbf{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^{-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 and $\tau$ 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}=\Lambda^{2}$. Here, $\Lambda$ is the real eigenvalue obtained using the a-FBPM simulations. Using this relation, we can obtain the lifetime of the phonon mode as

$$\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 consequently, $\tau$ can be sensitive to the mesh density. Hence, it is important to select an optimal grid to obtain the converged values of $\Lambda$. 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% when $N_{x}\geq 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:

$$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 an integral multiple of the time taken for each round trip, $t_{\text{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_{\text{YIG}}=200\mu$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 $\tau^{wophase}=0.1$ ms and $\tau^{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, $\tau^{wphase}=14.1$ ms, closely matches with $\tau$ 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: $\tau=\frac{-2t_{\text{HBAR}}}{v_{2}\text{ln}(E_{in}/E_{tot})}$. 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}\propto\int_{V_{in}}|\textbf{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 cavities [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.