Composition-tunable magnon-polaron anomalies in spin Seebeck effects
in epitaxial Bi_xY_3-xFe₅O₁₂ films

Figure 3 shows the cross-sectional TEM images of the Pt/YIG/GGG, Pt/Bi_0.5Y_2.5IG/GGG, and Pt/Bi_0.9Y_2.1IG/SGGG samples. The overall sample cross-section images [Figs. 3(a), 3(d), and 3(g)] show high uniformity and flatness of the films. The high resolution images at the Bi_xY_3-xIG/(S)GGG(111) interfaces [Figs. 3(c), 3(f), and 3(i)] reveal atomically smooth and epitaxial interfaces of our garnet films. Neither macroscopic defects nor misalignment in the lattice planes were observed in TEM images. As shown in the insets of Figs. 3(a), 3(d), and 3(g), the diffraction patterns of the Bi_xY_3-xIG films coincide with those of the (S)GGG substrates, indicating that the Bi_xY_3-xIG layers are grown as single-crystalline films with the [111] orientation in the out-of-plane direction. We also confirmed clear interfaces between the (polycrystalline) Pt-film and Bi_xY_3-xIG-film contacts [see Figs. 3(b), 3(e), and 3(h)], which allows us to investigate interfacial spin-current transport in these sample systems.

Figures 4(a), 4(b), and 4(c) display the $2\theta$-$\omega$ XRD patterns of the Bi_xY_3-xIG films around the (444) Bragg peaks from the GGG ($2\theta=51.058^{\circ}$) or SGGG ($2\theta=50.546^{\circ}$) substrates. For all the samples, we observed the (444) diffraction peaks from the Bi_xY_3-xIG layer, confirming the out-of-plane orientation of Bi_xY_3-xIG[111]/(S)GGG[111], consistent with the TEM results shown in Fig. 3. The XRD rocking curves exhibit full width at half maximum (FWHM) values of $0.0085^{\circ}$, $0.0093^{\circ}$, and $0.0103^{\circ}$ for the YIG, Bi_0.5Y_2.5IG, and Bi_0.9Y_2.1IG, respectively, which are close to the values for the substrates ($0.0078^{\circ}$ for GGG and $0.0099^{\circ}$ for SGGG), ensuring high crystalline quality of our films.

To further characterize structural properties of the films, we performed reciprocal space mapping (RSM) Kubota2012APEX ; Kubota2013JMMM ; Avci2017NatMat . Figures 4(d), 4(e), 4(f) show the RSM data around the (664) diffraction for the YIG/GGG, Bi_0.5Y_2.5IG/GGG, and Bi_0.9Y_2.1IG/SGGG samples, respectively. For all the samples, the asymmetric diffraction peaks lie on one vertical line, indicating that the films have a pseudomorphic structure with the in-plane lattice constant identical to that of the substrate along the $\langle 11\overline{2}\rangle$ direction Kubota2012APEX ; Kubota2013JMMM ; Avci2017NatMat . The pseudomorphic structure was confirmed also along the $\langle 1\overline{1}0\rangle$ direction by the RSM data for the (486) diffraction [see Figs. 4(g), 4(h), and 4(i)].

From the XRD and RSM results, the out-of-plane (in-plane) lattice constant was estimated to be $12.376~{}{\rm\AA}$ ($12.383~{}{\rm\AA}$) for the YIG film on GGG, $12.456~{}{\rm\AA}$ ($12.383~{}{\rm\AA}$) for the Bi_0.5Y_2.5IG film on GGG, and $12.446~{}{\rm\AA}$ ($12.508~{}{\rm\AA}$) for the Bi_0.9Y_2.1IG film on SGGG. This result shows that the YIG film grown on GGG is nearly free from strain, while the Bi_0.5Y_2.5IG (Bi_0.9Y_2.1IG) film grown on GGG (SGGG) has compressive (tensile) epitaxial strain. In TABLE 1, we show the in-plane biaxial strain $\epsilon_{||}$ and out-of-plane uniaxial strain $\epsilon_{\perp}$ values evaluated according to Ref. Kubota2013JMMM, . The in-plane stress $\sigma_{||}$ values Kubota2013JMMM are also calculated using the elastic coefficients $C_{11}\,(=\rho c_{\rm LA}^{2}$; $\rho$: mass density) and $C_{44}\,(=\rho c_{\rm TA}^{2})$ and the relation $C_{11}-C_{12}=2C_{44}$ Gurevich-Melkov_text ; YIG_Paoletti_text , which we list also in TABLE 1.

In Figs. 5(a), 5(b), and 5(c), we show the magnetic field $H$ dependence of the magnetization $M$ of the (a) YIG, (b) Bi_0.5Y_2.5IG, and (c) Bi_0.9Y_2.1IG films, respectively. Here, the $M$-$H$ curves are recorded at $300~{}\textrm{K}$ under the applied $H$ parallel to the in-plane $[11\overline{2}]$ (solid lines) and to the out-of-plane $[111]$ (dashed lines) axes. The pure YIG film (nearly free from strain) has an in-plane easy axis due to the shape anisotropy and negligibly-small magnetocrystalline anisotropy Kubota2013JMMM ; SSE_Kikkawa2013PRL ; Gallagher2016APL ; the saturation field $H_{\rm s}$ value for the in-plane $[11\overline{2}]$ field ($\mu_{0}H_{\rm s}^{||}=3~{}\textrm{mT}$) is much smaller than that for the out-of-plane $[111]$ field ($\mu_{0}H_{\rm s}^{\perp}=0.22~{}\textrm{T}$) [see Fig. 5(a)]. The in-plane saturation field $H_{\rm s}^{||}$ for the Bi_0.5Y_2.5IG (Bi_0.9Y_2.1IG) film is estimated to be $\mu_{0}H_{\rm s}^{||}=1~{}\textrm{mT}$ (0.17 T), which is much smaller (slightly higher) than that for the out-of-plane $\mu_{0}H_{\rm s}^{\perp}=0.23~{}\textrm{T}$ ($0.12~{}\textrm{T}$) [see Figs. 5(b) and 5(c)]. The different $H_{\rm s}^{||}$ v.s $H_{\rm s}^{\perp}$ features can be interpreted in terms of the stress-induced anisotropy; for a material with a negative magnetostriction constant $\lambda_{111}<0$ such as Bi_xY_3-xIG, the magnetic easy axis tends to lie in (perpendicular to) the film plane for compressive (tensile) epitaxial strain Hansen1983PRB ; Kubota2012APEX ; Kubota2013JMMM ; Tang2016PRB ; Fu2017APL . Therefore, for the Bi_0.5Y_2.5IG films on GGG with the compressive epitaxial strain, both the shape anisotropy and stress-induced anisotropy make the magnetic easy axis within an in-plane direction ($H_{\rm s}^{||}\ll H_{\rm s}^{\perp}$). In contrast, for the Bi_0.9Y_2.1IG film on SGGG with the tensile epitaxial strain, the stress-induced out-of-plane anisotropy appears and prevails against the shape anisotropy, which renders the magnetic easy axis perpendicular to the film plane ($H_{\rm s}^{||}>H_{\rm s}^{\perp}$) Soumah2018NatCommun .

The saturation magnetization $M_{\rm s}$ values were estimated to be $3.34~{}\mu_{\textrm{B}}$ ($0.164~{}\textrm{T}$), $3.62~{}\mu_{\textrm{B}}$ ($0.177~{}\textrm{T}$), and $3.66~{}\mu_{\textrm{B}}$ ($0.175~{}\textrm{T}$) at $300~{}\textrm{K}$, for the YIG, Bi_0.5Y_2.5IG, and Bi_0.9Y_2.1IG films, respectively, in units of the Bohr magneton (Tesla) (see Fig. 5). The $M_{\rm s}$ value increases by the Bi substitution, which is consistent with the previous reports Hansen1983PRB ; Hansen1984ThinSolidFilms . Figures 5(d), 5(e), and 5(f) show the $T$ dependence of $M_{\rm s}$ of the (d) YIG, (e) Bi_0.5Y_2.5IG, and (f) Bi_0.9Y_2.1IG films, respectively, for the $H$ direction parallel to the in-plane $[11\overline{2}]$ axis. For all the samples, the $M_{\rm s}$ value increases with decreasing $T$ and approaches to $\sim 5~{}\mu_{\rm B}$ at the lowest $T$, showing good agreement with the theoretical and previous experimental results YIG_Gilleo-Geller ; Hansen1983PRB ; Hansen1984ThinSolidFilms ; Kikkawa2017PRB .

To determine the sound velocities of the Bi_xY_3-xIG films, we employed a pump-and-probe magneto-optical imaging method combined with a Fourier transform (FT) process (see Refs. Hashimoto2017NatCommun, and Hashimoto2018PRB, for the details of the optical system in this experiment). When a pumped laser pulse with a spot radius of $\sim 2.5~{}\mu\textrm{m}$ is irradiated to the sample, elastic waves are excited and propagate radially from the excitation point. Via the magnetoelastic interaction, the elastic waves tilt the magnetic moments from an in-plane to out-of-plane direction that can be detected as a change of the Faraday rotation angle of the probe pulse Hashimoto2017NatCommun ; Hashimoto2018PRB . Figure 6(a) shows a snapshot of the temporal evolution of out-of-plane magnetization distribution of the Bi_0.9Y_2.1IG sample obtained at the time delay of $t=6.8$ ns between the pump and probe pulses. Two ring structures are visible; the inner (outer) ring is attributed to the coupled TA (LA) and spin waves having a smaller (larger) velocity. By Fourier transforming the observed propagating waveform with respect to the time and spatial coordinates, we obtain coupled phonon and spin waves dispersion relations. As shown in Fig. 6(b), two $k$-linear dispersions with different slopes are discerned; the dispersion showing the lower (larger) slope is assigned to the TA (LA) phonon branch, from which we can estimate the sound velocities of each phonon mode. Measurements were performed also for the Bi_0.5Y_2.5IG sample. We note that the crossings between phonon and spin wave branches are inaccessible, since their positions in the momentum direction are large and therefore out of range in the present experiment.

Significant decreases in the sound velocities for the Bi:YIG films are indeed confirmed through the above measurements and analysis. As shown in Table 1, the TA- and LA-phonon sound velocities of the Bi_0.5Y_2.5IG films were determined to be $c_{\rm TA}=3.34\times 10^{3}~{}\textrm{m/s}$ and $c_{\rm LA}=6.42\times 10^{3}~{}\textrm{m/s}$, giving values smaller than those of YIG. For the Bi_0.9Y_2.1IG film, $c_{\rm TA}=3.30\times 10^{3}~{}\textrm{m/s}$ and $c_{\rm LA}=6.28\times 10^{3}~{}\textrm{m/s}$, showing the smallest values among our samples.

Using the magnetostriction constant $\lambda_{111}$ values for $\mu$m-thick Bi_xY_3-xIG films [$\lambda_{111}=-2.819\cdot(1+0.75x)\times 10^{-6}$] shown in Ref. Soumah2018NatCommun, and the sound velocities $c_{\rm TA}\,(=\sqrt{C_{44}/\rho})$ in Table 1, the magnetoelastic coupling constant $B_{\perp}/2\pi\,[=-3\,a_{\rm Bi:YIG}^{3}\lambda_{111}C_{44}/(2\pi\hbar)\,]$ ($a_{\rm Bi:YIG}$: lattice constant) Kittel1949RevModPhys ; Gurevich-Melkov_text ; YIG_Paoletti_text was estimated to be $1.87$, $2.05$, and $2.57~{}\textrm{THz}$ for the YIG, Bi_0.5Y_2.5IG, and Bi_0.9Y_2.1IG films, respectively. The $B_{\perp}$ value monotonically increases by increasing the Bi amount $x$, which may be attributed to the enhanced spin-orbit coupling induced by the substitution of Bi³⁺ ions Kumar2019JPCM .

Now, we present the experimental results on the LSSE in the Pt/Bi_xY_3-xIG samples. Figure 7(a) shows the $H$ dependence of the LSSE coefficient $S$ of the Pt/YIG sample measured at $T=50~{}\textrm{K}$. In this conventional system, we observed a clear SSE signal as well as two magnon-polaron-induced spikes on top of the almost flat $S$ background signal [see the blue and red filled triangles in Fig. 7(a) and magnified view shown in Fig. 7(d)]. The peak at the lower (higher) $H$ appears at $\mu_{0}H_{\rm TA}=2.51~{}\textrm{T}$ ($\mu_{0}H_{\rm LA}=9.27~{}\textrm{T}$), consistent with Ref. Kikkawa2016PRL, . The $H_{\rm TA}$ and $H_{\rm LA}$ values coincide with the $H$ intensities at which the magnon dispersion curve touches the TA- and LA-phonon dispersion curves of YIG, respectively, being the conditions that the magnon and phonon modes can be coupled over the largest volume in momentum space, so the magnon-polaron effects are maximized Kikkawa2016PRL ; Flebus2017PRB .

Next, let us focus on the effect of Bi substitution on the LSSE and magnon-polaron features. Figures 7(b) and 7(c) show the $H$ dependence of $S$ of the Pt/Bi_0.5Y_2.5IG and Pt/Bi_0.9Y_2.1IG samples, respectively. We observed SSE signals in these samples, of which the amplitudes are comparable to that for the Pt/YIG shown in Fig. 7(a). The peak features show up also in these Pt/Bi:YIG systems as marked by the blue and red filled triangles in Figs. 7(b) and 7(c). Interestingly, the peak positions shift toward lower $H$ values by increasing the amount of Bi substitution. The peak positions at the lower and higher $H$ for the Pt/Bi_0.5Y_2.5IG sample are estimated to be $\mu_{0}H_{\rm TA}=2.15~{}\textrm{T}$ and $\mu_{0}H_{\rm LA}=8.18~{}\textrm{T}$, respectively. For the Pt/Bi_0.9Y_2.1IG, the peak fields are further decreased: $\mu_{0}H_{\rm TA}=2.06~{}\textrm{T}$ and $\mu_{0}H_{\rm LA}=7.47~{}\textrm{T}$.

The observed peak shifts are explained in terms of the reduction of the sound velocities and resultant touching fields of Bi:YIG films. The touching field can be evaluated quantitatively by solving the combined equation for the magnon ($\omega_{\rm mag}$) and phonon ($\omega_{\rm ph}$) dispersion relations and their group velocities:

where the magnon dispersion relation, disregarding the dipolar interaction ($M_{\rm s}\ll H_{\rm TA,LA}$), reads $\omega_{\rm mag}=D_{\rm ex}k^{2}+\gamma\mu_{0}H$, while the phonon dispersions are given by $\omega_{\rm ph}=c_{\rm TA,LA}k$. Here, $D_{\rm ex}$, $k$, and $\gamma$ represent the exchange stiffness, wavenumber, and gyromagnetic ratio, respectively. Calculation of Eq. (1) leads to

which shows that the touching field $H_{\rm TA(LA)}$ scales quadratically with the sound velocity $c_{\rm TA(LA)}$, and so does the anomaly field in the SSE. In fact, the observed peak fields $H_{\rm TA,LA}$ are well reproduced by the sound velocities $c_{\rm TA(LA)}$ shown in Table 1 and reasonable $D_{\rm ex}$ values of $8.0\times 10^{-6}$, $7.1\times 10^{-6}$, and $7.3\times 10^{-6}~{}\textrm{m}^{2}/\textrm{s}$ for the YIG, Bi_0.5Y_2.5IG, and Bi_0.9Y_2.1IG, respectively. The frequency value for the touching point $\omega_{\rm MTA}/2\pi$ ($\omega_{\rm MLA}/2\pi$) between the magnon and TA(LA)-phonon branches at $H_{\rm TA}$ ($H_{\rm LA}$), depicted as green filled circles in Fig. 1(b), can be evaluated as $\omega_{\rm MTA}/2\pi=0.15~{}\textrm{THz}$ ($\omega_{\rm MLA}/2\pi=0.52~{}\textrm{THz}$) for the YIG, $0.13~{}\textrm{THz}$ ($0.47~{}\textrm{THz}$) for the Bi_0.5Y_2.5IG, and $0.12~{}\textrm{THz}$ ($0.43~{}\textrm{THz}$) for the Bi_0.9Y_2.1IG.

We carried out systematic measurements of the temperature dependence of the LSSE in these Pt/Bi_xY_3-xIG samples. Figures 8(a), 8(d), and 8(g) show the $S$-$H$ curves for the various average sample temperature $T_{\mathrm{avg}}$ (from $300$ to $3~{}\textrm{K}$) for the Pt/YIG, Pt/Bi_0.5Y_2.5IG, and Pt/Bi_0.9Y_2.1IG samples, respectively. In all the samples, the overall $T$ dependences of $S$ agree with those for the Pt/YIG-film samples reported in Refs. Kikkawa2015PRB, and Jin2015PRB, ; when the sample temperature is decreased from $300~{}\textrm{K}$, the magnitude of $S$ increases and shows a broad peak at around $T=150~{}\textrm{K}$ [see Figs. 9(a)-9(c)]. On decreasing $T$ further, the $S$ signal begins to decrease and eventually goes to zero.

We now discuss the $T$ dependence of the magnon-polaron features in the LSSE. First, let us focus on the results at a low-$T$ range below $50~{}\textrm{K}$, at which characteristic $T$ response is expected to occur due to the competition between the thermal energy $k_{\rm B}T$ and the magnon-polaron’s touching energy $\hbar\omega_{\mathrm{MTA,MLA}}$, which governs the $T$ dependence of thermal occupation of magnon-polaron modes Kikkawa2016PRL . Figures 9(d), 9(e), and 9(f) represent the $T$ dependence of the magnon-polaron peak amplitude $\delta S_{\rm TA,LA}\equiv S(H_{\rm TA,LA})-S_{\rm TA0,LA0}$ for the Pt/YIG, Pt/Bi_0.5Y_2.5IG, and Pt/Bi_0.9Y_2.1IG samples, respectively, where $S_{\rm TA0(LA0)}$ is the extrapolated background $S$ at the peak position $H_{\rm TA(LA)}$ [see also the definition shown in Fig. 7(d) and its caption]. The peak intensities $\delta S_{\rm TA,LA}$ at $H=H_{\mathrm{TA}}$ and $H_{\mathrm{LA}}$ exhibit different $T$ dependences. We further introduce and hereafter discuss the quantity $\delta S_{\mathrm{TA(LA)}}/S_{\mathrm{TA0(LA0)}}$, where $\delta S_{\rm TA(LA)}$ is normalized by the background $S_{\rm TA0(LA0)}$, which allows us to evaluate the $T$-dependence of the magnon-polaron SSE coefficient relative to the background magnonic one at $H=H_{\mathrm{TA(LA)}}$ and also to exclude common $T$-dependent factors in the magnon-polaron and magnonic SSE voltages such as the interfacial spin conductance, electric resistance, and spin diffusion length in the Pt film Cornelissen2016PRB-chemical-potential ; Guo2016PRX ; Cornelissen2017PRB . As shown in Figs. 9(g)-9(i), $\delta S_{\mathrm{TA}}/S_{\mathrm{TA0}}$ monotonically increases with decreasing $T$ and takes a maximum value at the lowest $T$ for all the samples (see also the $S$-$H$ curves shown in Fig. 8). In contrast, $\delta S_{\mathrm{LA}}/S_{\mathrm{LA0}}$, is gradually suppressed below $T^{\ast}\sim 12~{}\text{K}$ ($\sim 7~{}\text{K}$) for the Pt/YIG sample (Pt/Bi_0.5Y_2.5IG and Pt/Bi_0.9Y_2.1IG samples) [see Figs. 8 and 9(g)-9(i)]. At the lowest $T$, the $S$ peak at $H_{\mathrm{LA}}$ is so small and almost indistinguishable from the background [see Figs. 8(c), 8(f), and 8(i)]. The $T$-dependent feature can be interpreted in terms of the difference in the frequency values of the branch touching point $\omega_{\mathrm{MTA}}$ for $H=H_{\mathrm{TA}}$ and $\omega_{\mathrm{MLA}}$ for $H=H_{\mathrm{LA}}$ [see Fig. 1(b)] Kikkawa2016PRL ; Flebus2017PRB . Here, in the unit of temperature, the $\omega_{\mathrm{MLA}}$ values approximately correspond to $T_{\rm MLA}(=\hbar\omega_{\rm MLA}/k_{\rm B})=25~{}\textrm{K}$, $22~{}\textrm{K}$, and $21~{}\textrm{K}$ for the YIG, Bi_0.5Y_2.5IG, and Bi_0.9Y_2.1IG, respectively, and they are more than 3 times greater than those of the $\omega_{\mathrm{MTA}}$ values ($7~{}\textrm{K}$, $6~{}\textrm{K}$, and $6~{}\textrm{K}$, respectively). Therefore, for $T<T_{\rm MLA}$, the excitation of magnon polarons at the touching frequency $\omega\sim\omega_{\mathrm{MLA}}$ is rapidly suppressed, which leads to the disappearance of the $S$ peak at $H_{\mathrm{LA}}$ at the lowest $T$. Furthermore, the reduction of $\omega_{\mathrm{MLA}}$ values for Bi_0.5Y_2.5IG and Bi_0.9Y_2.1IG, by a value of $\sim 5~{}\textrm{K}$, compared to that for the YIG may be responsible for the difference of the onset $T^{\ast}$ values between the Pt/Bi:YIG ($T^{\ast}\sim 7~{}\textrm{K}$) and Pt/YIG ($T^{\ast}\sim 12~{}\textrm{K}$), at which $\delta S_{\mathrm{LA}}/S_{\mathrm{LA0}}$ starts to be suppressed with decreasing $T$.

Next, we move on to the results at a higher-$T$ range ($50~{}\textrm{K}\lesssim T_{\rm avg}\lesssim 300~{}\textrm{K}$). In this $T$ range, we are also able to resolve small but finite peak structures at $H=H_{\mathrm{TA}}$ for all the samples, while, at $H_{\mathrm{LA}}$, the peak structures were detectable below $\sim 150~{}\textrm{K}$ ($\sim 250~{}\textrm{K}$) for the Pt/YIG (Pt/Bi:YIG) sample (see Figs. 8 and 9). As shown in Figs. 9(g)-9(i), the peak amplitudes at $H_{\mathrm{TA}}$ and $H_{\mathrm{LA}}$ relative to the background SSE signals, $\delta S_{\mathrm{TA}}/S_{\mathrm{TA0}}$ and $\delta S_{\mathrm{LA}}/S_{\mathrm{LA0}}$, monotonically decrease with increasing $T$, which may be due to the appearance of large contributions of thermally populated (uncoupled) magnons with the energy of $k_{\rm B}T>\hbar\omega_{\rm MTA,MLA}$ to the background SSE signal. Furthermore, the observed structures with only peak shapes (no dip shapes) suggest that the scattering rate of magnons $\tau^{-1}_{\rm mag}$ with the frequencies relevant to the hybridization (i.e., $\omega\sim\omega_{\rm MTA,MLA}$) are always larger than that of phonons $\tau^{-1}_{\rm ph}$ at all the $T$ range: $\tau^{-1}_{\rm mag}>\tau^{-1}_{\rm ph}$ Kikkawa2016PRL ; Flebus2017PRB . We also found that, in this $T$ range, as the temperature is increased, the $\delta S_{\mathrm{LA}}/S_{\mathrm{LA0}}$ value decreases much faster than the $\delta S_{\mathrm{TA}}/S_{\mathrm{TA0}}$ value for all the samples [see Figs. 9(g)-9(i)]. We note that with increasing $T$, dominant scattering sources for magnons and phonons may change from disorders (impurities) to inelastic magnon-magnon and phonon-phonon scatterings due to their increased thermal population Rezende2014PRB ; Cornelissen2016PRB-chemical-potential ; Boona2014PRB ; Shi2021PRL_YIG-bulk ; Schmidt2018PRB . A detailed knowledge on the $T$ and $k$ dependences of such scattering events for the Bi_xY_3-xIG films will thus be important to fully understand the observed behavior.

Finally, we compare the amplitudes of the peak structures in the LSSE in the Pt/Bi_xY_3-xIG samples. The magnon-polaron contribution $\delta S_{\mathrm{TA}}/S_{\mathrm{TA0}}$ at $H=H_{\rm TA}$ for the Pt/YIG, Pt/Bi_0.5Y_2.5IG, and Pt/Bi_0.9Y_2.1IG samples at the lowest $T\sim 3~{}\textrm{K}$ ($T=50~{}\textrm{K}$) was evaluated as 2.4, 1.3, and 5.2 (0.15, 0.11, and 0.25), respectively. This nonmonotonic dependence with respect to the Bi amount is not simply explained by the magnetoelastic coupling constant $B_{\perp}$, since it monotonically increases by increasing the Bi substitution (see TABLE 1). Although larger $B_{\perp}$ values should be beneficial for the magnon-polaron formation and resultant SSE anomalies, our results suggest that the scattering ratio parameterized by $\eta=\tau^{-1}_{\rm mag}/\tau^{-1}_{\rm ph}$ may play a rather important role for the comparison of the magnon-polaron peak intensities between each sample. Furthermore, in the LSSE results, we found that the magnon-polaron contribution at $H=H_{\rm TA}$ is greater than that at $H_{\rm LA}$, i.e., $\delta S_{\mathrm{TA}}/S_{\mathrm{TA0}}>\delta S_{\mathrm{LA}}/S_{\mathrm{LA0}}$ for all the $T$ range and all the samples [see Figs. 9(g)-9(i)]. This tendency is compared with that for the nlSSE results in the next section.

We now show the results on the nlSSE in the Pt/Bi_xY_3-xIG/Pt systems. Figures 10(a), 10(b), and 10(c) show the $H$ dependence of the nlSSE coefficient ${\tilde{S}}$ of the Pt/YIG/Pt, Pt/Bi_0.5Y_2.5IG/Pt, and Pt/Bi_0.9Y_2.1IG/Pt samples measured at $T=50~{}\textrm{K}$, respectively. The overall signal appears with a sign opposite to that for the LSSE, indicating that thermally-excited magnons accumulate beneath the detector Pt strip. The observed sign is indeed a characteristic in nlSSEs for $\mu$m-thick garnet films with a long injector-detector separation distance $d$ ($d=8~{}\mu\textrm{m}$ in the present case) Cornelissen2017PRB ; Shan2016PRB ; Shan2017PRB . As marked by the blue and red unfilled triangles in the nlSSE signals in Fig. 10, magnon-polaron induced anomalies are also observed at the touching fields $H_{\rm TA,LA}$, which shift toward lower $H$ values by increasing the Bi amount [see also the magnified views for $H>0$ shown in Figs. 10(d), 10(e), and 10(f)]. The $H_{\rm TA(LA)}$ values are evaluated as $\mu_{0}H_{\rm TA}=2.48~{}\textrm{T}$ ($\mu_{0}H_{\rm LA}=9.17~{}\textrm{T}$) for the Pt/YIG/Pt, $2.13~{}\textrm{T}$ ($8.14~{}\textrm{T}$) for the Pt/Bi_0.5Y_2.5IG/Pt, and $2.02~{}\textrm{T}$ ($7.38~{}\textrm{T}$) for the Pt/Bi_0.9Y_2.1IG/Pt, in good agreement with those estimated from the LSSE results. For the nlSSE, the magnon-polaron formation decreases rather than increases the ${\tilde{S}}$ signal at the touching fields $H_{\rm TA,LA}$, causing dip structures (note that the background nlSSE signal is opposite in sign compared to that for the LSSECornelissen2017PRB ; Shan2018APL ; Oyanagi2020AIPAdv ). In Ref. Cornelissen2017PRB, , the dip feature is discussed in terms of the competition between a temperature-gradient induced magnon current, ${\bf J}_{\rm m}^{T}=-(\zeta/T)\nabla T$, from the injector to detector and a diffusive backflow, ${\bf J}_{\rm m}^{\mu}=-\sigma_{\rm m}\nabla\mu_{\rm m}$, induced by the gradient of the magnon chemical potential $\mu_{\rm m}$ [${\bf J}_{\rm m}^{T}\;||+{\bf\hat{z}}$ and ${\bf J}_{\rm m}^{\mu}\;||-{\bf\hat{z}}$ in the coordinate system shown in Fig. 2(b)], which are both increased by the hybridization with phonons when $\tau^{-1}_{\rm mag}>\tau^{-1}_{\rm ph}$ Cornelissen2017PRB ; Flebus2017PRB . In particular, assuming that the magnon spin conductivity $\sigma_{\rm m}$ is more strongly increased than the bulk SSE coefficient $\zeta$ via the magnon-polaron formation, the increased backflow of magnons toward the injector direction causes the reduction of magnon accumulation beneath the Pt detector, which leads to the suppression of the nlSSE at the touching fields Cornelissen2017PRB .

Through a careful look at the magnon-polaron anomalies in Fig. 10, we notice that the intensity of the anomaly at $H=H_{\rm LA}$ is larger than that at $H_{\rm TA}$ in the nlSSE for all the samples: $|\delta{\tilde{S}}_{\mathrm{TA}}/{\tilde{S}}_{\mathrm{TA0}}|<|\delta{\tilde{S}}_{\mathrm{LA}}/{\tilde{S}}_{\mathrm{LA0}}|$. The $\delta{\tilde{S}}_{\rm TA(LA)}/{\tilde{S}}_{\rm TA0(LA0)}$ values at $H_{\rm TA(LA)}$ are evaluated to be $-0.029$ ($-0.035$) for the Pt/YIG/Pt, $-0.014$ ($-0.077$) for the Pt/Bi_0.5Y_2.5IG/Pt, and $-0.079$ ($-0.16$) for the Pt/Bi_0.9Y_2.1IG/Pt at $T=50~{}\textrm{K}$, which we summarize in Fig. 11. The observed relationship for the nlSSE, $|\delta{\tilde{S}}_{\mathrm{TA}}/{\tilde{S}}_{\mathrm{TA0}}|<|\delta{\tilde{S}}_{\mathrm{LA}}/{\tilde{S}}_{\mathrm{LA0}}|$, holds in a wide $T$ range for $T\lesssim 100~{}\textrm{K}$, and is opposite to that for the LSSE results, where $\delta S_{\mathrm{TA}}/S_{\mathrm{TA0}}>\delta S_{\mathrm{LA}}/S_{\mathrm{LA0}}$ (see Fig. 11 and also Figs. S1 and S2 in the Supplemental Material SM for the $T$ dependence of the nlSSE). In the following, we discuss possible origins of the anisotropic feature. One possible scenario is a spectral non-uniform magnon current $J_{\rm m}(\omega)$ Kikkawa2015PRB ; Schmidt2021PRB that may vary depending on the experimental configurations; the frequency-resolved current $J_{\rm m}(\omega_{\rm MTA})$ at $H_{\rm TA}$ may be larger (smaller) than $J_{\rm m}(\omega_{\rm MLA})$ at $H_{\rm LA}$ for the longitudinal (nonlocal) configuration. In Refs. Kikkawa2015PRB, , Oyanagi2020AIPAdv, and Jin2015PRB, , through high-$H$ dependence experiments and its comparison with theory, a spectral non-uniformity in a magnon current is suggested to be present in both LSSEs and nlSSEs due to the frequency-dependent magnon scattering time Schmidt2021PRB ; Streib2019PRB . Moreover, in the present case, as shown in Figs. 7 and 10, the high-$H$ response of the background SSEs in the longitudinal and nonlocal configurations differs with each other. This may be a signature that the frequency-dependent $J_{\rm m}(\omega)$ contribution is different between the LSSE and nlSSE. Another possibility would be the anisotropy in magnon and magnon-polaron dispersion relations due to the nature of the dipole and magnetoelastic interaction Shen2015PRL ; Flebus2017PRB , which depend on the magnon propagation direction ${\bf J}_{\rm m}$ ($||~{}{\bf k}~{}||~{}\nabla T$) relative to the external magnetic field ${\bf H}$. In the LSSE, ${\bf J}_{\rm m}$ is essentially perpendicular to ${\bf H}$, while in the nlSSE parallel to ${\bf H}$ in the one-dimensional limit (see Fig. 2). Therefore, the LSSE and nlSSE may be affected by the anisotropic dispersion relations in a different manner. In fact, a solution of a Boltzmann transport equation shows that, for $\nabla T\perp{\bf H}$ (compatible with the LSSE), the magnon-polaron anomaly at $H=H_{\rm TA}$ is larger than that at $H_{\rm LA}$, while, for $\nabla T~{}||~{}{\bf H}$ (compatible with the nlSSE in the one-dimensional limit), the anomaly at $H_{\rm TA}$ is smaller than that at $H_{\rm LA}$ Flebus2017PRB . This is consistent with our experimental results, and suggests that the anisotropic nature of magnon and magnon-polaron dispersion relations may affect the SSEs. It is worth mentioning that the profile of the temperature gradient created by the local Pt heater in the nonlocal configuration may be affected by size effects Cornelissen2017PRB ; Shan2016PRB ; Shan2017PRB , such as the thickness of the magnetic layer and the size of the local heater. Future systematic measurements of the Bi:YIG thickness, heater size, and injector-detector separation distance ($d$) dependencies may provide useful information to understand possible roles of the size effects in the anisotropic feature of magnon-polaron signals. Besides, the magnon-polaron formation can also affect the magnon and phonon thermal conductivities Flebus2017PRB , which may modify the temperature gradient and resultant spin-current intensity at the onset field of magnon-polaron formation. The process, which has not been considered in analysis so far, would be important for further quantitative argument on the magnon-polaron anomalies in SSEs.