THz Control of Exchange Mode in a Ferrimagnetic Cavity

Abstract

The interaction of terahertz (THz) radiation with high-frequency spin resonances in complex magnetic materials is central for modern ultrafast magnonics. Here we demonstrate strong variations of the excitation efficiency of the sub-THz exchange magnon in a single crystal ferrimagnet (Gd,Bi)₃Fe₅O₁₂. An enhancement of the exchange magnon amplitude is observed when its frequency matches an eigenmode of the cavity created by the sample interfaces. Moreover, this enhancement is accompanied by a 5-fold decrease in effective damping of the exchange mode. The THz-exchange magnon interaction in the cavity is analyzed within the developed Landau-Lifshitz-Gilbert formalism for three coupled magnetization sublattices and cavity-enhanced THz field. This work presents a novel approach for the THz excitation of spin dynamics in ferrimagnets and outlines promising pathways for the controlled optimization of light-spin coupling in single crystals.

I INTRODUCTION

In an era where large language models are becoming increasingly dominant in society, the need for faster computational power is more pronounced than ever before [1]. However, the predominant form of data storage technology today is still ferromagnetic-based, which limits operational frequencies to a few GHz [2]. To achieve significantly higher frequencies, i.e. in the THz range, antiferromagnetic magnetic recording and control is the subject of intense investigation nowadays  [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. Among media with antiferromagnetic ordering, ferrimagnets are especially appealing as they combine the feasibility of magnetic control with high THz-scale eigenfrequencies. One material class at the center of these investigations has been the rare earth iron garnets [11, 12, 13, 15]. Iron garnets crystallize in the l$\bar{\text{a}}$3d structure where the iron ions form two sublattices sitting in oxygen octahedral $\vec{M}_{\text{a}}$ and tetrahedral sites $\vec{M}_{\text{d}}$ with a ratio of 10 $\mu_{\text{B}}$:15 $\mu_{\text{B}}$, while the rare earth ions occupy the interstitial sites. By substituting the main rare earth element [16], and by doping with non-magnetic ions [17], the magnetic properties of these materials become highly tunable. Due to the strong Fe-Fe exchange interaction, iron garnets typically are treated as ferrimagnets characterized by two sublattices: Gadolinium $\vec{M}_{\text{Gd}}=\vec{M}_{\text{c}}$ and joint iron $\vec{M}_{\text{Fe}}=\vec{M}_{\text{d}}+\vec{M}_{\text{a}}$, respectively. The strength of the exchange interaction between these sublattices can be directly probed by measuring the corresponding magnetic exchange mode [18]. The frequency of this resonance is highly sensitive to temperature, especially in the vicinity of the magnetic compensation point, where $\lvert\vec{M}_{\text{Fe}}\rvert$ = $\lvert\vec{M}_{\text{Gd}}\rvert$. Furthermore, substitution of Gd with Bi enhances the magneto-optical response, enabling ultrafast magneto-optical spectroscopy [19, 20]. The exchange mode can be excited either resonatly [21] or via opto-magnetic effects [11].

There are several approaches to increase the coupling strength between photons and magnons: some are based on energy concentration using cavities [22, 23, 24], while others exploit multiple terahertz pulses as a stimulus to coherently control the magnon response [25, 26]. In this paper, we combine the concepts of cavity engineering with multiple THz pulse excitation to demonstrate the vast tuning possibilities of the coupling strength between light and magnetization at THz frequencies. Specifically, we analyze the practically important case where the cavity is created by the two interfaces of a plane-parallel micron-scale plate made of the ferrimagnet itself. The tunability of the exchange mode frequency through temperature control allows us to couple it with multiple cavity modes in a controllable way. By performing THz-pump optical-probe experiments, we monitor the amplitude and the frequency of the exchange mode to reveal the corresponding magnon-photon coupling strength at various temperatures. The enhancement of the mode amplitude is readily observable when the frequency of the magnon mode is resonant with an eigenmode of the cavity, while a reduction of the effective damping is observed across the full temperature range between 80 and 130 K. Numerical simulation based on Landau-Lifschitz-Gilbert equations in a two-sublattice system in a cavity show good agreement with the experimental findings. Overall, this work touches on a previously unexplored realm of THz cavity spin dynamics and provides a systematic understanding of the dynamical magnetic response of a ferrimagnetic cavity.

II MATERIALS AND METHODS

A 200 $\mu$m thick single crystal of (Gd,Bi)₃Fe₅O₁₂ (GdBIG) oriented along the (111) direction was grown using liquid phase epitaxy and polished to the prescribed thickness [11]. The sample was mounted in a continuous flow cryostat and cooled using liquid nitrogen with a dc magnetic field applied in the sample plane. Broadband THz pulses were generated using the tilted-pulse-front optical rectification technique in a LiNbO₃ crystal [27]. For that optical pulses from a Ti:sapphire amplifier with 6 mJ per pulse, 1 kHz repetition rate and a center wavelength of 800 nm were used. Part off this beam was split of to be used as a low intensity probe pulse. Two types of THz spectroscopy experiments were conducted: THz time-domain (TDS) spectroscopy and THz-pump optical-probe (TPOP). In TDS, the THz beam was transmitted through the GdBIG sample, after which the THz pulse field strength was detected via electro-optic sampling in a 1 mm thick ZnTe crystal [28]. TPOP was performed by overlapping the 800 nm probe beam with the THz beam within the sample. The THz-induced changes of the probe polarization were detected using the balanced photodiode technique [28]. Both the THz and optical radiation were at normal incidence. For both experimental configurations, the polarization state and fluence of the THz radiation was controlled using two wire grid polarizers. In TDS the THz electric field strength was reduced below 10 $\frac{\text{kV}}{\text{cm}}$, whereas for TPOP, THz field strengths up to 700 $\frac{\text{kV}}{\text{cm}}$ were used.

III RESULTS

Figure 1: Static characterization of GdBIG. (a) THz electric field transmitted through the sample $E_{\text{Sam}}$ and an empty sample holder $E_{\text{Ref}}$ at room temperature, offset for better visibility. (b) Transmission spectrum of GdBIG at different temperatures. The reference spectrum is shown as a gray background, where the sharp absorption lines result from the water vapor absorption. The cavity modes are numbered according to their order.

As a first step, the transmission properties of GdBIG were studied by THz TDS spectroscopy at different temperatures. Figure 1(a) shows the THz electric field transmitted through the sample $E_{\text{Sam}}$ and the reference THz electric field $E_{\text{Ref}}$ at room temperature. The main peak of $E_{\text{Sam}}$ is delayed by $\sim$3.1 ps compared to the reference signal. Furthermore, a second peak at $\sim$10.7 ps is visible, originating in the back reflection of the THz pulse from the rear interface of GdBIG. For a sample thickness of 200 $\mu$m and the time delay observed in $E_{\text{Ref}}$ the effective refractive index $n$, in the $\sim 1$ THz frequency range, is estimated to be $\sim$5.6. The amplitude transmission spectra plotted in Fig. 1(b) are calculated as the ratio of the absolute values of the $E_{\text{Sam}}$ and $E_{\text{Sam}}$ FFT spectra. A strong periodic modulation of the transmission is observed, which is caused by formation of a THz pulse train due to multiple reflections from both crystal interfaces [29]. The peaks shift to lower frequencies with decreasing temperature indicating an increase in the refractive index. Moreover, the sample becomes non-transparent above 1.3 THz at 300 K and above 1.1 THz at 5 K, indicating the presence of an absorption line outside of the spectral window which is softening with decreasing temperature. This absorption line is absent in pure Gd₃Fe₅O₁₂ [30], which suggests that this could be a phonon mode involving Bi in the interstitial sites. This notion is reinforced by the appearance of a Raman-active optical phonon at 1.8 THz when substituting Gd with Tb in Gd_2.34Tb_0.66Fe₅O₁₂ [31]. Without changing the coupling strength, the mass ratio between the Tb and Bi ions would lower this phonon frequency to $\approx$1.6 THz. Furthermore, room temperature Raman measurements (Fig. S1) shows the presence of a peak at 1.28 THz. According to the mutual exclusion principle [32], the appearance of this mode in both Raman and IR spectra can be explained by the lowering of the crystal symmetry due to the random distribution of Bi dopants in GdBIG [33].

Figure 2: Temperature dependence of the Gd-Fe exchange mode. (a) THz-induced polarization rotation for 80 K, 100 K, and 130 K, offset for better visibility. (b) Fourier spectra of the time traces from panel (a), the arrows indicate the exchange mode, offset for better visibility. (c) The frequency (black circles) and the effective damping (red triangles) of the exchange mode as a function of temperature. Both values are extracted with a Lorentzian function fit to the exchange mode peak. The black line shows a fit based on the Kaplan-Kittel exchange resonance formula (see Eq. S5) [18] with the molecular field parameter $\text{N}_{\text{Gd-Fe}}=-0.80\frac{\text{mol}}{\text{cm}^{3}}$. The red line is a fit with $\alpha^{\text{eff}}_{\text{exch}}=\alpha^{\text{eff}}\left(|\vec{M}_{\text{Gd}}|+|\vec{M}_{\text{Fe}}|\right)/\left(|\vec{M}_{\text{Gd}}|-|\vec{M}_{\text{Fe}}|\right)$ [34] with the constant effective Gilbert damping factor of the sublattices $\alpha^{\text{eff}}=0.0046$.

Next, TPOP spectroscopy was performed. Figure 2(a) shows THz induced polarization rotation at three distinct temperatures of 130 K, 100 K, and 80 K. The time traces indicate a rich frequency spectrum which is strongly temperature-dependent. Corresponding Fourier spectra are plotted in Fig. 2(b) showing multiple peaks composing the induced polarization rotation in Fig. 2(a). Two of the most pronounced features, centered at 0.13 THz and 0.26 THz, do not move with temperature, while the frequency of the 0.4 THz peak at 130 K increases to 0.7 THz upon lowering the temperature down to 80 K. The extracted peak position with a temperature step of 2 K is shown in Fig. 2(c), in an excellent agreement with the temperature dependence observed for the exchange resonance between the Gd and Fe sublattices [35]. The positions of the other visible peaks in Fig. 2(b) that do not change with temperature within our measurement accuracy, are consistent with the 1st and the 2nd cavity modes in Fig. 1(b). Therefore, they can convincingly be attributed to the forced response of the spin system to the THz pulse, shaped by multiple reflections at GdBIG surfaces. Moreover, the cavity modes modify the observable line-width of the magnetic resonance, thereby enabling control of the apparent effective damping of the exchange mode, $\alpha^{\text{eff}}_{\text{exch}}$. The apparent change in damping results from the recurrent parametric excitation of the spin system by a phase-matched sequence of the THz pulses within the cavity. The extracted line-width of the exchange mode at every temperature is shown in Fig. 2(c). The value of $\alpha^{\text{eff}}_{\text{exch}}$ varies from 0.007 to 0.04 in the observed temperature range. This is almost one order of magnitude lower than previously reported measurements on similar compounds. With a non-resonant excitation the intrinsic damping of the exchange mode was shown to be from 0.05 to 0.1 [36]. Rigorous estimation of the enhancement is challenging due to limitations in frequency resolution and the low finesse of the cavity. However, we observed that the enhancement is consistent with the finesse, measuring 5-6 in our experiment.
Controlling the exchange mode necessitates understanding the underlying excitation mechanism. To this end, the THz field strength dependence of the induced magnetization dynamics was analyzed. As illustrated in Fig. 3(a), the amplitude scales linearly with the strength of the THz field. Moreover, the change of exchange mode amplitude with changing THz field polarization $\beta$, where $\beta=\angle{(\vec{M},\vec{H}_{\text{THz}})}$, shown in Fig. 3(b), reveals a sine-like dependence. The two observed maxima at $\beta=90^{\circ}$ and $\beta=-90^{\circ}$ are in agreement with the linear Zeeman torque excitation mechanism [25].

Figure 3: THz light - exchange mode coupling mechanism. The exchange mode amplitude as a function of THz pulse field strength (a) and THz field polarization (b). These measurements were performed at a temperature of 100 K. The solid lines are the expected dependencies for a Zeeman excitation mechanism. (c) Amplitude of the exchange mode as function of its frequency. The spectrum of the transmitted pulse train corresponding to the magnetic field of the THz pulse $H_{\text{THz}}$ is shown as a gray background. (d) Amplitude of the exchange mode as function of its frequency derived from the simulations (see Sec. IV). The simulated excitation spectrum of $H_{\text{THz}}$ is shown as a gray background.

Taking advantage of the high tunability of the exchange mode with temperature, we analyzed its coupling with the cavity modes and their influence on the resulting spin dynamics. Figure 3(c) shows the amplitude of the exchange mode as a function of its frequency, where each data point corresponding to a distinct temperature. Although the net magnetic moment of GdBIG significantly decreases as the temperature approaches the compensation point $T_{\text{comp}}=223$ K [11], this monotonic trend does not account for the observed periodic behavior. However, consistent with the Zeeman torque excitation mechanism, the amplitude of the exchange mode linearly depends on the THz field strength. To illustrate this, we present the spectrum of the magnetic field of the transmitted THz pulse, $H_{\text{THz}}\propto E_{\text{Sam}}$, measured with EOS at 100 K, shown as a gray background in Fig. 3(c). Indeed, the exchange mode amplitude increases twofold when its frequency matches one of the cavity resonances.

IV THEORETICAL ANALYSIS AND DISCUSSION

To understand the THz-induced magnetization dynamics in the presence of cavity modes, we employ the modified Landau–Lifshitz–Gilbert equation formalism. Taking into account the dominant role of the iron-iron exchange, the magnetization dynamics can be described by coupled equations for a two-sublattice ferrimagnet [37]:

$$\frac{d\vec{M}_{\text{i}}}{dt} =-(\mu_{0}\gamma_{\text{i}})\left[\vec{M}_{\text{i}}\times\vec{H}^{\text{eff}}_{\text{i}}\right]-(\mu_{0}\gamma_{\text{i}})\frac{\alpha_{\text{i}}}{\lvert\vec{M}_{\text{i}}\rvert}\vec{M}_{\text{i}}\times\left[\vec{M}_{\text{i}}\times\vec{H}^{\text{eff}}_{\text{i}}\right]$$ (1)

with $\mu_{0}$ the vacuum magnetic permeability, $M_{\text{i}}$ the net magnetizations, $\gamma_{\text{i}}$ the gyromagnetic ratios, and $\alpha_{\text{i}}$ the Gilbert damping factors of the Gd and Fe sublattice dynamics. $\vec{H}^{\text{eff}}_{\text{i}}$ is the effective magnetic field defined as

$$\vec{H}^{\text{eff}}_{\text{Gd,Fe}}=\vec{H}_{\text{app}}+\vec{H}_{\text{THz}}-N_{\text{Gd-Fe}}\vec{M}_{\text{Fe,Gd}}.$$ (2)

Here, the first term with $\mu_{0}\vec{H}_{\text{app}}=150$ mT represents the applied external magnetic field, $\vec{H}_{\text{THz}}$ denotes the THz pulse magnetic field, and $N_{\text{Gd-Fe}}$ the effective Gd-Fe exchange parameter. To model the impact of the cavity, we introduce multiple delayed THz pulses so that the temporal profile of the THz stimulus closely follows the experimentally observed transmitted field $E_{\text{Sam}}$ in Fig. 1(b) at 100 K. A comparison of the simulated and measured THz waveforms and spectra is presented in Fig. S2(a,b) of the Appendix B.

The exchange field which is responsible for the coupling of the two sublattices, is represented by the last term in Eq. 2. It is worth noting that the two-sublattice model fails to account for the temperature dependence of the effective Gd-Fe coupling. Instead, a three-sublattice model should be considered where the two Fe magnetizations are treated separately (see Appendix C). Moreover, by applying the three-sublattice model, the connection between $N_{\text{Gd-Fe}}$ and parent molecular field coefficients can be derived. From fitting the temperature dependence of the exchange mode frequency with the Kaplan-Kittel model [18], we obtain $N_{\text{Gd-Fe}}$ = -0.89 $\frac{\text{mol}}{\text{cm}^{3}}$, which is in a good agreement with the value derived from the SQUID data [35]. Away from the compensation point $T_{\text{comp}}=223$ K it is safe to assume that the Gilbert damping factors $\alpha_{\text{Fe}}=\alpha_{\text{Gd}}=\alpha$ [34] and almost temperature independent [36], while their values were set to 0.03 based on the previous measurements on this sample [11] and $\gamma_{\text{Fe}}\approx\gamma_{\text{Gd}}=\gamma$ = $28\frac{\text{GHz}}{\text{T}}$ [37].

By numerically solving Eqs. (1) using the Euler method, the magnetization dynamics of the sublattices was extracted. To compare the simulation results with the experiment it is important to note that the Faraday rotation at 800 nm, in our experimental geometry, is mostly sensitive to the magnetization of the iron sublattice in GdBIG [11]. Figure S2(c) shows the simulated and experimentally observed THz-induced magnetization dynamics at 100 K. Both time scans show a fast oscillation attributed to the exchange mode dynamics. The experimental scan has an additional low frequency component that is not captured by our simulations. However, the simulation does not take into account propagation effects and dispersion of the sample, therefore an exact match is not expected. Simulating temperature changes by varying the sublattice magnetizations allows us to extract the amplitude and the frequency of the exchange mode shown in Fig. 3(d). The gray background illustrates the spectrum of the applied THz magnetic field $H_{\text{THz}}$. Similar to the experiment, the amplitude of the exchange mode is enhanced when its frequency coincides with the cavity modes created by the sample interfaces.

We believe that the concept of driving high-frequency spin dynamics with a THz pulse train in a cavity formed by the crystal interfaces can be further explored. For example, we suggest a sandwich structure consisting of a THz generation crystal and a ferrimagnet capped by optically transparent and highly THz reflective interfaces. The progress in THz coating makes it possible to achieve very high reflection coefficients, which accentuate the action of the THz pulse train on the medium. For example, indium–tin–oxide-coated glass has THz amplitude reflection coefficient $r$ as high as 0.9 [38], while using THz Bragg mirrors this value potentially can reach 0.99 [39]. This is substantially higher than that obtained in a bare GdBIG crystal with Fresnel reflection at interfaces, where we estimate $r$ = 0.7. To study the impact of interface reflection, we simulate THz-induced magnetization dynamics with $r$ = 0.95 and plot the amplitude as a function of the frequency of the exchange mode in Fig. S3(a). A much higher dynamical range of the mode amplitude variations is observed, while $\alpha^{\text{eff}}_{\text{exch}}$ can be further reduced by fourfold. This is illustrated in Fig. S3(b) by comparing magnetization dynamics waveforms at a temperature of 128 K. Moreover, in conjunction with using ferrimagnets with much lower intrinsic damping, e.g., yttrium iron garnet [40], exchange mode coupling with THz light might even enter the strongly coupled regime  [41, 42].

V CONCLUSION

In conclusion, we investigated the coupling between light at THz frequencies and the Gd-Fe spin exchange mode in GdBIG using THz-pump optical-probe measurements. We show that tuning the exchange mode into a resonance with the cavity modes strongly enhances the efficiency of the THz excitation of high-frequency spin dynamics. The excitation mechanism is attributed to the Zeeman torque. Moreover, the cavity reduces the effective damping of the exchange mode 5-fould. The results are analyzed by means of numerical simulations based on a system of coupled Landau-Lifshitz-Gilbert equations which show a good agreement with the experimental observations. A three-subsystem approach developed here for the first time allowed to conclude that the effective Gd-Fe molecular field strength exhibits no variations in the 70-130 K temperature range. This work improves the understanding of efficient control of the exchange modes at THz frequencies. Furthermore, it opens an intuitively clear and uncomplicated path to achieve strong coupling between light and sub-THz exchange modes by means of designing the effictive interfaces of a magnetic material.

The supporting data and codes for this article are available from Zenodo [43].

VI Appendix

VI.1 Raman Measurement

The Raman spectrum was recorded by a triple-grating spectrometer using a 532 nm continuous wave laser in backscattering geometry at a power of 1 mW on the sample. The direction of polarisation of the analysed light was parallel to that of the incoming beam. The Raman spectrum in the relevant spectral range is shown in Fig. S1.

Figure S1: Raman intensity measured on GdBIG in parallel polarization configuration at room temperature.

VI.2 Modeling of THz Magnetic Field

Figure S2: (a) Magnetic THz field applied in the simulation compared to the experimentally measured THz field transmitted through the sample at 100 K. (b) Corresponding Fourier spectra. (c) Comparison of the THz-induced magnetization dynamics in the sample between experiment and simulation at 100 K.

Figure S3: (a) Numerical simulation of the behavior of the exchange mode amplitude as a function of its frequency with two different reflection coefficients at the interface. (b) THz-induced magnetization dynamics modeled with the reflection coefficient $r=0.7$ (bottom curve) and $r=0.95$ (top curve).

To model the THz pulse magnetic field $H_{\text{THz}}$, we used sequences of pulses defined as

$$H_{\text{THz}}(t)=\sum_{k=0}^{N}H_{0}\cdot r^{2k}\text{sin}(\omega_{0}t-k\cdot\omega_{0}\tau)\cdot e^{-(t-k\cdot\tau)^{2}/2\sigma^{2}},$$ (S1)

with $H_{0}$ the magnetic field amplitude, $\omega_{0}$ the center frequency, $\sigma$ the pulse duration, $N$ the number of reflections, $r$ amplitude reflection coefficient at the crystal interface, and $\tau$ the time delay introduced by propagation in the medium. To follow the shape of the THz pulse used in the experiment we choose $\mu_{0}H_{0}=0.266$ T, $\omega_{0}/2\pi=0.6$ THz, and $\sigma=0.63$ ps. The coefficients $r$ and $\tau$ are calculated using a constant refractive index of $n=5.6$ extracted from the THz-TDS measurements in Fig. 1. $N=10$ is selected to ensure that all reflections within the chosen time window are accounted for. Figure S2(a) shows the resulting sequences of THz pulses compared to the THz electric field measured with EOS after transmitting through the sample at 100 K. A very good agreement of the corresponding spectra is also evident, see Fig. S2(b). To compare the THz-induced magnetization dynamics with that observed in the experiment we plot the out-of-plane component of $\vec{M}_{\text{Fe}}$ in Fig. S2(c).

This approach allows for the testing of different paths to optimize the excitation efficiency of the exchange mode. Figure S3 shows the simulation results of THz-induced spin dynamics with the reflection coefficients at the interfaces of $r$ = 0.95. The other parameters are kept the same. It is seen that improving the cavity finesse up to 60 results in a threefold enhancement of the dynamic range of the spin exchange mode peak amplitude (green circles).

VI.3 Effective exchange constant

Rare earth iron garnets are typically considered two-sublattice magnets. This assumption is based on the fact that the molecular field coefficient N${}_{\text{ad}}$, which accounts for the interaction between iron at the octahedral and tetrahedral sites, is much larger than the iron-rare earth (rare earth is indexed by ”c”) molecular field coefficients N${}_{\text{cd}}$ and N${}_{\text{ca}}$ [44]. For that reason, the majority of recent papers do not address the connection between the effective coupling constant N${}_{\text{Gd-Fe}}$ and the three-sublattice molecular field coefficients. However, this assumption potentially overlooks the temperature dependence of the effective coupling constant [45]. Considering a three-sublattice model instead provides an additional way to estimate N${}_{\text{Gd-Fe}}$ by analysing SQUID data in combination with molecular field theory. The main steps for deriving N${}_{\text{Gd-Fe}}$ are provided below.

To address the dynamics of a three-sublattice magnet, we write three Landau–Lifshitz equations assuming no damping for simplicity:

$$\frac{d\vec{M}_{\text{a}}}{dt}=-(\mu_{0}\gamma_{\text{Fe}})\left[\vec{M}_{\text{a}}\times\left(\vec{H}_{\text{THz}}-N_{\text{ad}}\vec{M}_{\text{d}}-N_{\text{ac}}\vec{M}_{\text{c}}\right)\right]$$

$$\frac{d\vec{M}_{\text{d}}}{dt}=-(\mu_{0}\gamma_{\text{Fe}})\left[\vec{M}_{\text{d}}\times\left(\vec{H}_{\text{THz}}-N_{\text{da}}\vec{M}_{\text{a}}-N_{\text{dc}}\vec{M}_{\text{c}}\right)\right]$$ (S2)

$$\frac{d\vec{M}_{\text{c}}}{dt}=-(\mu_{0}\gamma_{\text{Gd}})\left[\vec{M}_{\text{c}}\times\left(\vec{H}_{\text{THz}}-N_{\text{ca}}\vec{M}_{\text{a}}-N_{\text{cd}}\vec{M}_{\text{d}}\right)\right]$$

Then, we look for a solution in the form $\vec{M}_{i}=\vec{M}_{i}^{\text{init}}+\vec{m}_{i}e^{i2\pi ft}$, where $\vec{M}_{i}^{\text{init}}$ are the ground state magnetizations, $\vec{m}_{i}$ are the corresponding perturbations, and $i$ corresponds to the indices $a$, $d$ and $c$. We assume that all magnetizations are collinearly aligned in a ground state, note that $N_{\text{ad}}=N_{\text{da}}$, $N_{\text{ca}}=N_{\text{ac}}$, and $N_{\text{cd}}=N_{\text{dc}}$, take into account terms only linear on $\vec{m}_{i}$ and exclude terms of the second order on N${}_{\text{ac}}$ and N${}_{\text{dc}}$. Then, the Gd-Fe exchange frequency can be found:

$$f_{\text{Gd-Fe}}=\frac{\mu_{0}\gamma}{2}\left[N_{\text{ad}}\left(M_{d}-M_{a}\right)-N_{\text{ac}}\left(M_{a}+M_{c}\right)+N_{\text{dc}}\left(M_{d}-M_{c}\right)\right]-\\ \frac{\mu_{0}\gamma}{2}\left[N_{\text{ad}}^{2}\left(M_{a}^{2}+M_{d}^{2}-2M_{a}M_{d}\right)-2N_{\text{ac}}N_{\text{ad}}\left(M_{a}^{2}+M_{a}M_{c}-M_{a}M_{d}+M_{c}M_{d}\right)-2N_{\text{ad}}N_{\text{dc}}\left(M_{d}^{2}-M_{a}M_{d}-M_{a}M_{c}-M_{c}M_{d}\right)\right]^{1/2}$$ (S3)

$$f_{\text{Gd-Fe}}=\frac{\mu_{0}\gamma}{2}\left[N_{\text{ad}}\left(M_{d}-M_{a}\right)-N_{\text{ac}}\left(M_{a}+M_{c}\right)+N_{\text{dc}}\left(M_{d}-M_{c}\right)\right]-\\ \frac{\mu_{0}\gamma}{2}N_{\text{ad}}\left(M_{d}-M_{a}\right)\left[1-2\left(\frac{N_{\text{ac}}}{N_{\text{ad}}}\right)\frac{M_{a}^{2}+M_{a}M_{c}-M_{a}M_{d}+M_{c}M_{d}}{M_{a}^{2}+M_{d}^{2}-2M_{a}M_{d}}-2\left(\frac{N_{\text{dc}}}{N_{\text{ad}}}\right)\frac{M_{d}^{2}-M_{a}M_{d}-M_{a}M_{c}-M_{c}M_{d}}{M_{a}^{2}+M_{d}^{2}-2M_{a}M_{d}}\right]^{1/2}$$ (S4)

Given that the iron-iron exchange interaction considerably exceeds the rare earth-iron exchange interaction, $N_{\text{ac}}$, $N_{\text{dc}}<<N_{\text{ad}}$, the last two terms inside the square root in Eq. S4 represent a small parameter to which a Taylor expansion can be applied. This results in

$$f_{\text{Gd-Fe}}=-(\mu_{0}\gamma)N_{\text{Gd-Fe}}\left|M_{\text{d}}-M_{\text{a}}-M_{\text{Gd}}\right|;$$ (S5)

$$N_{\text{Gd-Fe}}=-{\frac{N_{\text{cd}}M_{\text{d}}+N_{\text{ca}}M_{\text{a}}}{M_{\text{d}}-M_{\text{a}}}}.$$ (S6)

Equation S6 clearly shows that since the magnetizations $\vec{M}_{\text{d}}$ and $\vec{M}_{\text{a}}$ are temperature-dependent, the effective Gd-Fe exchange coupling also exhibits temperature dependence when the three-sublattice model is considered. Substituting molecular field coefficients N${}_{\text{cd}}$ = 6 $\frac{\text{mol}}{\text{cm}^{3}}$ and N${}_{\text{ca}}$ = -3.44 $\frac{\text{mol}}{\text{cm}^{3}}$ from [46] gives estimation of N${}_{\text{Gd-Fe}}$ changing in the range from -0.89 to -0.90 $\frac{\text{mol}}{\text{cm}^{3}}$ in the studied temperature range from 80 to 130 K.