Acoustically driven magnon-phonon coupling in a layered antiferromagnet

First, IDTs (35 nm aluminium) and electrodes (5 nm titanium / 200 nm gold) are deposited onto 128$\degree$ Y-cut LiNbO₃ chips. The IDT fingers are 400 nm wide with 1.2 $\upmu$m spacing, giving a SAW wavelength 3.2 $\upmu$m and frequency 1.1 GHz. The distance between IDT1 and IDT2 is approximately 600 $\upmu$m. Next, bulk CrCl₃ is exfoliated onto polydimethylsiloxane (PDMS) sheets using sticky tape (Nitto). Flakes with uniform thickness are transferred onto LiNbO₃ between IDT1 and IDT2 using a conventional PDMS dry stamping technique. Bulk CrCl₃ crystals are obtained from the commercial suppliers 2D Semiconductors (USA) and HQ Graphene (Netherlands).

The LiNbO₃ chip is mounted on a radio-frequency compatible chip carrier and loaded into either a Montana closed-cycle cryostat with external electromagnet in 1 axis (Sample 1), or a helium bath cryostat with superconducting magnet coils in 2 axes (Sample 2). The former has a base temperature around 5 K and the latter around 4.2 K. Both cryostats allow variable sample temperature up to at least 30 K. Coaxial cables are used to connect the chip carrier to a vector network analyzer which is capable of measuring SAW transmission at 1.1 GHz. A time gating function is applied to the signal in order to filter out electromagnetic noise and retrieve the signals S21 and S12 at longer timescales 150 - 250 ns.

We consider a smiple Heisenberg model of antiferromagnet

$$H=-J_{\parallel}\sum_{\langle m,n\rangle}\left(\bm{S}^{A}_{m}\cdot\bm{S}^{A}_{n}+\bm{S}_{m}^{B}\cdot\bm{S}^{B}_{n}\right)-J_{\perp}\sum_{\{m,n\}}\bm{S}^{A}_{m}\cdot\bm{S}_{n}^{B},$$

(4)

where $\bm{S}^{A}_{n},\bm{S}^{B}_{n}$ denote spins belonging to $A$ and $B$ sublattices respectively with $n$ labelling the lattice sites. For modeling CrCl₃, we take $J_{\parallel}>0$ the ferromagnetic intra-layer exchange and $J_{\perp}<0$ the antiferromagnetic inter-layer exchange interactions respectively while $\langle m,n\rangle$ and $\{m,n\}$ denote intra- and inter-sublattice nearest neighbour links.

We use the mean-field ansatz

$$\langle\bm{S}_{n}^{A}\rangle=\langle S\rangle\hat{\bm{z}},\quad\langle\bm{S}_{n}^{B}\rangle=-\langle S\rangle\hat{\bm{z}},$$

(5)

which gives the molecular fields

$$\bm{B}_{A}=\left(z_{\parallel}J_{\parallel}-z_{\perp}J_{\perp}\right)\langle S\rangle\hat{\bm{z}},\quad\bm{B}_{B}=-\left(z_{\parallel}J_{\parallel}-z_{\perp}J_{\perp}\right)\langle S\rangle\hat{\bm{z}},$$

(6)

where $z_{\parallel}$ and $z_{\perp}$ are the numbers of intra- and inter-sublattice links per atom. The expectation value $\langle S\rangle$ is determined by solving

$$\frac{\langle S\rangle}{S}=B_{S}\left(q\frac{\langle S\rangle}{S}\right),\quad q=\frac{z_{\parallel}J_{\parallel}-z_{\perp}J_{\perp}}{k_{B}T},$$

(7)

where $B_{S}\left(x\right)$ is the Brillouin function

$$B_{S}\left(x\right)=\frac{2S+1}{2S}\coth\frac{2S+1}{2S}x-\frac{1}{2S}\coth\frac{1}{2S}x.$$

(8)

The asymptotic expansion

$$B_{S}\left(x\right)=\frac{S+1}{3S}x-\frac{S+1}{3S}\frac{2S^{2}+2S+1}{30S^{2}}x^{3}+\cdots,\quad x\rightarrow 0$$

(9)

implies at $T=T_{N}$

$$q_{T_{N}}=\frac{3S}{S+1}=\frac{z_{\parallel}J_{\parallel}-z_{\perp}J_{\perp}}{k_{B}T_{N}}$$

(10)

so that one can eliminate the microscopic parameters in favour of $T_{N}$;

$$q=\frac{3S}{S+1}\frac{T_{N}}{T}.$$

(11)

Throughout all the calculations in this work, we take $T_{N}=14$ K. For the spin parameter $S$, we have two choices; the nominal spin of Cr³⁺ $S=3/2$, or the macrospin approximation $S\rightarrow\infty$. Since past literature report a two-step phase transition where the 2D honeycomb layers first order ferromagnetically, which is followed by the antiferromagnetic order in the out-of-plane direction McGuire et al. (2017), we think the latter is more appropriate and use it to compute

$$M_{s}\left(T\right)=M_{s}\left(0\right)\lim_{S\rightarrow\infty}\frac{\langle S\rangle}{S},$$

(12)

where $\langle S\rangle/S$ on the right-hand-side is taken to be the self-consistent (numerical) solution of Eq. (7). We note that the other choice $S=3/2$ gives a similar temperature dependence, and the difference is irrelevant considering the qualitative nature of our analysis. In this formulation, $\bm{B}_{A}=-\bm{B}_{B}$ should be proportional to the exchange field $H_{E}$ appearing in the spin wave analysis, which yields

$$H_{E}\left(T\right)=H_{E}\left(0\right)\lim_{S\rightarrow\infty}\frac{\langle S\rangle}{S}.$$

(13)

We fixed the values of $M_{s}\left(0\right)$ and $H_{E}\left(0\right)$ for Sample 1 such that $M_{s}\left(T=1.56~{}{\rm K}\right)=250$ mT and $H_{E}\left(T=1.56~{}{\rm K}\right)=100$ mT and used them to generate Fig. 1d,e in the main text. Note that the reference temperature of 1.56 K was chosen so as to facilitate comparison with Ref. Macneill et al. (2019).

For Sample 2, the uniaxial anisotropy also appears to be important. We model it by adding the following term to the macroscopic free energy density;

$$F_{u}=-K_{u}\left\{\left(\hat{\bm{u}}\cdot\bm{n}_{A}\right)^{2}+\left(\hat{\bm{u}}\cdot\bm{n}_{B}\right)^{2}\right\},$$

(14)

where $\hat{\bm{u}}$ is the unit vector along the easy axis, and $K_{u}$ is the strength of anisotropy in units of energy density (for definitions of free energy density and $\bm{n}_{A,B}$, see the next section). Defined in this phenomenological way, the temperature dependence of $K_{u}$ is well established theoretically in a seminal work by Callen Callen and Callen (1966) to be $K_{u}\propto M_{s}\left(T\right)^{3}$. In computing the transmission spectra, this form was assumed and resulted in the temperature dependence of Fig. 3.

Although the focus of the present work is magneto-elastic coupling, a large part of the experimental results can be understood by considering only magnetic properties of CrCl₃. Let $\bm{n}_{A}=\bm{M}_{A}/M_{s}$ and $\bm{n}_{B}=\bm{M}_{B}/M_{s}$ be the normalized magnetization vectors for the respective sublattice, and introduce spherical coordinate variables by

$$\bm{n}_{A}=\begin{pmatrix}\sin\theta_{A}\cos\phi_{A}\\ \sin\theta_{A}\sin\phi_{A}\\ \cos\theta_{A}\\ \end{pmatrix},\quad\bm{n}_{B}=\begin{pmatrix}\sin\theta_{B}\cos\phi_{B}\\ \sin\theta_{B}\sin\phi_{B}\\ \cos\theta_{B}\\ \end{pmatrix}.$$

(15)

We set our coordinate $z$-axis to be along the crystal $c$-axis of CrCl₃. All the macroscopic magnetic properties should be derivable from an appropriately constructed free energy density $F$. For our purposes, it is sufficient to assume the following form:

	$\displaystyle F$	$\displaystyle=$	$\displaystyle J_{E}\left\{\sin\theta_{A}\sin\theta_{B}\cos\left(\phi_{A}-\phi_{B}\right)+\cos\theta_{A}\cos\theta_{B}\right\}-K_{\perp}\left(\cos^{2}\theta_{A}+\cos^{2}\theta_{B}\right)$		(16)
			$\displaystyle-\frac{K_{u}}{2}\left\{\sin^{2}\theta_{A}\cos 2\left(\phi_{A}-\phi_{u}\right)+\sin^{2}\theta_{B}\cos 2\left(\phi_{B}-\phi_{u}\right)\right\}-\mu_{0}M_{s}\bm{H}\cdot\left(\bm{n}_{A}+\bm{n}_{B}\right).$

Here $J_{E}>0$ is the antiferromagnetic exchange energy density (not to be confused with the microscopic exchange energies $J_{\parallel},J_{\perp}$ in the previous section), $K_{\perp}$ is the out-of-plane uniaxial anisotropy that arises from the intra-layer demagnetizing field and spin-orbit interactions, $K_{u}\geq 0$ is an externally induced in-plane uniaxial anisotropy that breaks the 6-fold rotation symmetry of CrCl₃, $\phi_{u}$ represents its easy-axis direction that is taken to be a free parameter, and $\bm{H}$ is the external magnetic field. While we assume that the out-of-plane anisotropy is dominated by the demagnetizing field $K_{\perp}=-M_{s}^{2}/2$, the spin-orbit contribution might not be entirely negligible. Although including it can change the theoretical temperature dependence, corrections to the mean-field approximation is far more likely sources of discrepancy so that we do not pursue this direction any further.

We take $\bm{H}$ to be in the $ab$-plane and parameterize it by

$$\bm{H}=H\begin{pmatrix}\cos\phi\\ -\sin\phi\\ 0\\ \end{pmatrix},\quad H>0.$$

(17)

The unusual sign convention for the $y$-component is in accordance with the clockwise convention of the polar plots in the main text. The equilibrium magnetization configuration is determined by minimizing $F$. It is clear that $\theta_{A}=\theta_{B}=\pi/2$. While we speak of spin waves, the wavelength relevant to our study is of order 1 $\mu$m, for which the effect of exchange interactions is expected to be subdominant. Therefore, we treat them as if they were spatially uniform modes. The linearized Landau-Lifshitz equation reads

$$\left\{i\frac{\omega}{\gamma\mu_{0}}\begin{pmatrix}0&-1&0&0\\ 1&0&0&0\\ 0&0&0&-1\\ 0&0&1&0\\ \end{pmatrix}-\begin{pmatrix}A_{1}&0&C_{1}&0\\ 0&A_{2}&0&C_{2}\\ C_{1}&0&B_{1}&0\\ 0&C_{2}&0&B_{2}\\ \end{pmatrix}\right\}\begin{pmatrix}\delta\theta_{A}\\ \delta\phi_{A}\sin\theta_{A}\\ \delta\theta_{B}\\ \delta\phi_{B}\sin\theta_{B}\\ \end{pmatrix}=0,$$

(18)

where $\theta_{A},\phi_{A},\theta_{B},\phi_{B}$ now refer to the equilibrium state and $\delta\theta_{A},\delta\phi_{A},\delta\theta_{B},\delta\phi_{B}$ small perturbations around it, and

	$\displaystyle A_{1}$	$\displaystyle=$	$\displaystyle-H_{E}\cos\left(\phi_{A}-\phi_{B}\right)+H\cos\left(\phi+\phi_{A}\right)+M_{s}+H_{u}\cos 2\left(\phi_{A}-\phi_{u}\right),$
	$\displaystyle A_{2}$	$\displaystyle=$	$\displaystyle-H_{E}\cos\left(\phi_{A}-\phi_{B}\right)+H\cos\left(\phi+\phi_{A}\right)+2H_{u}\cos 2\left(\phi_{A}-\phi_{u}\right),$
	$\displaystyle B_{1}$	$\displaystyle=$	$\displaystyle-H_{E}\cos\left(\phi_{A}-\phi_{B}\right)+H\cos\left(\phi+\phi_{B}\right)+M_{s}+H_{u}\cos 2\left(\phi_{B}-\phi_{u}\right),$
	$\displaystyle B_{2}$	$\displaystyle=$	$\displaystyle-H_{E}\cos\left(\phi_{A}-\phi_{B}\right)+H\cos\left(\phi+\phi_{B}\right)+2H_{u}\cos 2\left(\phi_{B}-\phi_{u}\right),$
	$\displaystyle C_{1}$	$\displaystyle=$	$\displaystyle H_{E},$
	$\displaystyle C_{2}$	$\displaystyle=$	$\displaystyle H_{E}\cos\left(\phi_{A}-\phi_{B}\right).$

The eigenfrequencies are obtained to be

	$\displaystyle\frac{\omega^{2}}{\gamma^{2}\mu_{0}^{2}}$	$\displaystyle=$	$\displaystyle\frac{A_{1}A_{2}+B_{1}B_{2}+2C_{1}C_{2}}{2}\pm\Bigg{\{}\left(\frac{A_{1}A_{2}-B_{1}B_{2}}{2}\right)^{2}+\frac{A_{1}B_{1}+A_{2}B_{2}}{2}\left(C_{1}^{2}+C_{2}^{2}\right)$		(19)
			$\displaystyle+\left(A_{1}A_{2}+B_{1}B_{2}\right)C_{1}C_{2}-\frac{A_{1}B_{1}-A_{2}B_{2}}{2}\left(C_{1}^{2}-C_{2}^{2}\right)\Bigg{\}}^{1/2}$
		$\displaystyle=$	$\displaystyle\frac{A_{1}A_{2}+B_{1}B_{2}+2C_{1}C_{2}}{2}$
			$\displaystyle\pm\sqrt{\left(\frac{A_{1}A_{2}-B_{1}B_{2}}{2}\right)^{2}+A_{1}B_{1}C_{2}^{2}+A_{2}B_{2}C_{1}^{2}+\left(A_{1}A_{2}+B_{1}B_{2}\right)C_{1}C_{2}}.$

Setting $H_{u}=0$, they reduce to the frequency equivalent of Eq. (1) in the main text. In generating Fig. 3a in the main text, we minimized $F$ in Eq. (16) numerically to obtain $\phi_{A},\phi_{B},\theta_{A},\theta_{B}$, and evaluated Eq. (19) with the minus sign to plot the acoustic mode resonance frequencies. The parameters for Sample 2, which are temperature dependant via Eqs. (12) and (13), were set at $T=1.56$ K as $\mu_{0}M_{s}=250$ mT, $\mu_{0}H_{E}=130$ mT, $\mu_{0}H_{u}=2.4$ mT, and $\phi_{u}=-\pi/20$ rad $\approx-9\degree$. Note that $H_{E}$ being different from Sample 1 is not surprising considering that the two samples were taken from crystals grown in different conditions.

In this subsection, we discuss the coupling between Rayleigh surface acoustic wave and the antiferromagnetic spin waves described in the previous section. Let an isotropic elastic body (the substrate plus the magnetic film on top) occupy the half space $z<0$ and assume there is no stress applied on the surface $z=0$. Acoustic waves in an isotropic media are characterized by only two parameters; longitudinal and transverse sound velocities $c_{P}$ and $c_{S}$. When Rayleigh surface acoustic wave with frequency $\omega=c_{R}k$ is propagating along $x$-axis, the displacement vector $\bm{u}$ is given by

$$\bm{u}=\Re\left[C\begin{pmatrix}\left(1-\xi_{S}^{2}\right)\left\{e^{\kappa_{P}z}-\left(1-\xi_{S}^{2}\right)e^{\kappa_{S}z}\right\}\\ 0\\ -i\sqrt{1-\xi_{P}^{2}}\left\{\left(1-\xi_{S}^{2}\right)e^{\kappa_{P}z}-e^{\kappa_{S}z}\right\}\\ \end{pmatrix}e^{-i\left(\omega t-kx\right)}\right].$$

(20)

Here $c_{R}$ is the Rayleigh wave velocity solely determined by $c_{P},c_{S}$, $C$ is a constant, and

$$\kappa_{P}=k\sqrt{1-\frac{c_{R}^{2}}{c_{P}^{2}}},\quad\kappa_{S}=k\sqrt{1-\frac{c_{R}^{2}}{c_{S}^{2}}},\quad\xi_{P}^{2}=\frac{c_{R}^{2}}{c_{P}^{2}},\quad\xi_{S}^{2}=\frac{c_{R}^{2}}{2c_{S}^{2}}.$$

(21)

The physical discussions should be based on the strain tensor $\epsilon_{ab}=\left(\partial_{b}u_{a}+\partial_{a}u_{b}\right)/2$ instead of $\bm{u}$ itself. The boundary condition enforces $\epsilon_{zx}=0$ at the surface and it stays close to zero within $\sim 1/k$ from the surface. In our setup, for the $\sim 100$ nm films, we can therefore assume $\epsilon_{zx}$ is subdominant. Then the only nonzero components of the strain tensor to be taken into account are $\epsilon_{xx}$ and $\epsilon_{zz}$.

We introduce the free energy density of magneto-elasticity to derive magnon-phonon coupling. Assuming full rotational symmetry, we use the following form:

$$F_{\rm me}=b\epsilon_{ab}\left(n_{a}^{A}n_{b}^{A}+n_{a}^{B}n_{b}^{B}\right)+2c\epsilon_{ab}n_{a}^{A}n_{b}^{B}.$$

(22)

$b$ corresponds to the usual magneto-elastic coupling of ferromagnetic materials while the inter-sublattice coefficient $c$ is peculiar to antiferromagnetic materials. Before going into the detailed calculation, let us see what kind of angular dependence one should expect for the magnon-SAW coupling. First of all, we are interested in the linear dynamics for which the free energy should be quadratic in the small fluctuations. The strain tensor $\epsilon_{ab}$ itself is already a small fluctuation, so that we need to keep only the first order terms in the perturbation of the magnetizations. Denoting the perturbations $\delta n_{a}^{A},\delta n_{a}^{B}$ and using $n_{a}^{A},n_{b}^{B}$ to refer to the fixed ground state values, the quadratic free energy reads

$$F_{\rm me}\approx 2b\epsilon_{ab}\left(n_{a}^{A}\delta n_{b}^{A}+n_{a}^{B}\delta n_{b}^{B}\right)+2c\epsilon_{ab}\left(n_{a}^{A}\delta n_{b}^{B}+n_{a}^{B}\delta n_{b}^{A}\right).$$

Next, as far as Rayleigh surface acoustic waves coupled to a thin magnetic field are concerned, as discussed above, we will need to keep only $\epsilon_{xx}$ and $\epsilon_{zz}$. However, because we consider only the ground states in the plane, $n_{z}^{A}=n_{z}^{B}=0$. Therefore, one can reduce the free energy further to obtain

$$F_{\rm me}\sim 2\epsilon_{xx}\left\{b\left(n_{x}^{A}\delta n_{x}^{A}+n_{x}^{B}\delta n_{x}^{B}\right)+c\left(n_{x}^{A}\delta n_{x}^{B}+n_{x}^{B}\delta n_{x}^{A}\right)\right\}.$$

(23)

Let us emphasize that $n_{x}^{A},n_{x}^{B}$ are not dynamical variables in this expression but just fixed coefficients that depend on $H$ and $\phi$, which is the root cause of angle dependence of the magnon-phonon couplings. The situation is slightly more complicated, however, since $\bm{n}^{A}$ and $\delta\bm{n}^{A}$ (and similarly for $B$) are orthogonal to each other so that some extra angle dependence arises from $\delta n_{x}^{A},\delta n_{x}^{B}$. To be quantitative, we write the perturbed magnetisation vectors as

	$\displaystyle\bm{n}^{A}$	$\displaystyle=$	$\displaystyle\begin{pmatrix}\cos\phi_{A}\\ \sin\phi_{A}\\ 0\\ \end{pmatrix},\quad\delta\bm{n}^{A}\ =\ \delta\phi_{A}\begin{pmatrix}-\sin\phi_{A}\\ \cos\phi_{A}\\ 0\\ \end{pmatrix}-\delta\theta_{A}\begin{pmatrix}0\\ 0\\ 1\\ \end{pmatrix},$
	$\displaystyle\bm{n}^{B}$	$\displaystyle=$	$\displaystyle\begin{pmatrix}\cos\phi_{B}\\ \sin\phi_{B}\\ 0\\ \end{pmatrix},\quad\delta\bm{n}^{B}\ =\ \delta\phi_{B}\begin{pmatrix}-\sin\phi_{B}\\ \cos\phi_{B}\\ 0\\ \end{pmatrix}-\delta\theta_{B}\begin{pmatrix}0\\ 0\\ 1\\ \end{pmatrix}.$

Note that the ground state angles $\phi_{A},\phi_{B}$ appear not only in the ground state direction but also multiplying the in-plane fluctuations $\delta\phi_{A},\delta\phi_{B}$. Substituting these into Eq. (23) yields

	$\displaystyle F_{\rm me}$	$\displaystyle\sim$	$\displaystyle-2\epsilon_{xx}\Big{\{}b\left(\delta\phi_{A}\cos\phi_{A}\sin\phi_{A}+\delta\phi_{B}\cos\phi_{B}\sin\phi_{B}\right)$		(24)
			$\displaystyle+c\left(\delta\phi_{A}\cos\phi_{B}\sin\phi_{A}+\delta\phi_{B}\cos\phi_{A}\sin\phi_{B}\right)\Big{\}}.$

If it were a ferromagnet, the direction of magnetization would roughly follow the magnetic field $\phi_{A}\approx\phi$ so that this expression explains why the magnon-phonon coupling is proportional to $\sin 2\phi=2\cos\phi\sin\phi$ and maximised at $\phi=45\degree$. Since we are dealing with an antiferromagnet, $\phi_{A},\phi_{B}$ have more complicated relations with $\phi$. In addition, the eigenmodes are acoustic and optical (only approximately if $H_{u}\neq 0$), i.e. in-phase and out-of-phase precessions of $\delta\bm{n}^{A}$ and $\delta\bm{n}^{B}$. Thus let us introduce new variables

$$\delta\phi_{\rm ac}=\delta\phi_{A}+\delta\phi_{B},\quad\delta\phi_{\rm op}=\delta\phi_{A}-\delta\phi_{B}.$$

(25)

In terms of those eigenmode variables, Eq. (24) reads

	$\displaystyle F_{\rm me}$	$\displaystyle\sim$	$\displaystyle-\epsilon_{xx}\Big{[}b\left\{\delta\phi_{\rm ac}\sin\left(\phi_{A}+\phi_{B}\right)\cos\left(\phi_{A}-\phi_{B}\right)+\delta\phi_{\rm op}\sin\left(\phi_{A}-\phi_{B}\right)\cos\left(\phi_{A}+\phi_{B}\right)\right\}$		(26)
			$\displaystyle+c\left\{\delta\phi_{\rm ac}\sin\left(\phi_{A}+\phi_{B}\right)+\delta\phi_{\rm op}\sin\left(\phi_{A}-\phi_{B}\right)\right\}\Big{]},$

where we used some trigonometric identities to simplify the result. Therefore, in order to understand the angular dependence of the magnon-phonon coupling in antiferromagnets, one needs to know $\phi_{A}+\phi_{B}$ and $\phi_{A}-\phi_{B}$ as a function of $\phi$. They are in general complicated. However, if there is no in-plane anisotropy, by symmetry considerations, we expect (note our convention of $\phi$ in Eq. (17))

$$\phi_{A}+\phi=-\left(\phi_{B}+\phi\right)\equiv\frac{\pi}{2}-\phi_{\rm cant}.$$

(27)

The notation is based on the following observation: In the limit of weak field $H\rightarrow 0$, $\bm{n}^{A}$ and $\bm{n}^{B}$ are antiparallel to each other and perpendicular to $\bm{H}$ so that the canting angle $\phi_{\rm cant}=0$. $\phi_{\rm cant}$ should monotonically increase as $H$ increases. With this, one obtains

$$\phi_{A}+\phi_{B}=-2\phi,\quad\phi_{A}-\phi_{B}=\pi-2\phi_{\rm cant}\left(H\right).$$

(28)

Importantly, $\phi_{A}-\phi_{B}$ is independent of $\phi$. Thus, one derives

	$\displaystyle g_{\rm ac}$	$\displaystyle\propto$	$\displaystyle-\left\{b\cos\left(\phi_{A}-\phi_{B}\right)+c\right\}\sin\left(\phi_{A}+\phi_{B}\right)\ =\ -\left(b\cos 2\phi_{\rm cant}-c\right)\sin 2\phi,$		(29)
	$\displaystyle g_{\rm op}$	$\displaystyle\propto$	$\displaystyle-\left\{b\cos\left(\phi_{A}+\phi_{B}\right)+c\right\}\sin\left(\phi_{A}-\phi_{B}\right)\ =\ -\left(b\cos 2\phi+c\right)\sin 2\phi_{\rm cant}.$		(30)

Therefore, the acoustic mode coupling is proportional to $\sin 2\phi$ while the optical coupling is $\cos 2\phi$ for the intra-sublattice term $\propto b$ and angle independent for the inter-sublattice term $\propto c$.

In order to calculate the SAW transmission amplitude, one needs to be more systematic. Following the approach taken by Refs. Verba et al. (2019); Yamamoto et al. (2022) for magnon-SAW coupling in ferromagnetic materials, one may reduce the equations of motion to the following form:

	$\displaystyle\left\{i\frac{\omega}{\gamma\mu_{0}}\begin{pmatrix}0&-1&0&0\\ 1&0&0&0\\ 0&0&0&-1\\ 0&0&1&0\\ \end{pmatrix}-\begin{pmatrix}A_{1}&0&C_{1}&0\\ 0&A_{2}&0&C_{2}\\ C_{1}&0&B_{1}&0\\ 0&C_{2}&0&B_{2}\\ \end{pmatrix}\right\}\begin{pmatrix}\delta\theta_{A}\\ \delta\phi_{A}\sin\theta_{A}\\ \delta\theta_{B}\\ \delta\phi_{B}\sin\theta_{B}\\ \end{pmatrix}=\epsilon_{R}\begin{pmatrix}0\\ g_{A}\\ 0\\ g_{A}\\ \end{pmatrix},$		(31)
	$\displaystyle\left\{\rho\left(\frac{\omega^{2}}{k^{2}}-c_{R}^{2}\right)+i\eta\omega\right\}\epsilon_{R}=\sigma+\overline{g_{A}}\delta\phi_{A}\sin\theta_{A}+\overline{g_{B}}\delta\phi_{B}\sin\theta_{B},$		(32)

where $\epsilon_{R}$ is an appropriately normalised amplitude of SAW, $\rho=4650$ kg/m³ is the mass density of LiNbO₃, $c_{R}\sim 3800$ m/s is the velocity of Rayleigh mode, $\sigma$ is the external stress generated by the input IDT, and $\eta$ is the coefficient of viscosity in the Kelvin-Voight model of viscoelasticity Rose (1999). $g_{A}$ and $g_{B}$ are the effective magnon-phonon coupling coefficients in the thin film limit arising from the isotropic magneto-elastic interaction (22):

	$\displaystyle g_{A}$	$\displaystyle=$	$\displaystyle-iC_{R}\sqrt{kd}\left[b\sin 2\phi_{A}+c\left\{\sin\left(\phi_{A}-\phi_{B}\right)+\sin\left(\phi_{A}+\phi_{B}\right)\right\}\right],$		(33)
	$\displaystyle g_{B}$	$\displaystyle=$	$\displaystyle-iC_{R}\sqrt{kd}\left[b\sin 2\phi_{B}+c\left\{\sin\left(\phi_{B}-\phi_{A}\right)+\sin\left(\phi_{A}+\phi_{B}\right)\right\}\right],$		(34)

where $C_{R}$ is a constant of order unity Yamamoto2020. We analytically solved Eqs. (31) and (32) for $\epsilon_{R}$ with $\omega=c_{R}k$, $b=0,C_{R}\sqrt{kd}c=10^{6}$, numerically computed $\phi_{A},\phi_{B}$ for given $H,\phi$, and evaluated $\epsilon_{R}$ to generate Figs. 3 and 4.

As a closing remark, we note that the theoretical model here is meant for capturing qualitative trends. In particular, the theory predicts very sharp lines for optical modes, while the experimental data point to multiple peak structure with a large broadening. There are two main factors that may cause the disagreement:

1. 1.

   The relative height of acoustic and optical peaks depends strongly on the precise form of magneto-elastic free energy, which may contain many more terms than included in Eq. (22.

2. 2.

   The broadening arising from fluctuations and inhomogeneity is not accounted for, which can become important when the optical mode frequency nears zero.

The growth quality of commercially obtained van der Waals magnetic materials is currently quite poor, but in time the material quality will likely improve, bringing with it our understanding of the above points. For the sake of presentation, the color coding in the simulated polar plots of Figs. 3 and 4 is based on a biased normalization. We use dB units (logarithmic scale), and assign the brightest color to a value of transmission appropriately large compared with the actual minimum of the data set. This is appropriate given that the figures intend to display the magnetic resonance positions, rather than amplitudes.