I n 1973, it was proposed 1,2 that photon and matter fields spontaneously appear in thermal equilibrium as a static transverse electromagnetic field and a static polarization, respectively, when the photon-matter coupling strength exceeds a certain threshold, entering the so-called ultrastrong coupling regime [3][4][5] . This phenomenon is known as the superradiant phase transition (SRPT) or Dicke phase transition, as the Dicke model was used in the theoretical calculations 1,2 , having been originally developed to describe the superradiance phenomena 6 .

The realization of the SRPT in thermal equilibrium may be expected to provide a new avenue for decoherence-robust quantum technology because the ground state of the Dicke model provides a quantum-squeezed vacuum on a photon-atom two-mode basis [7][8][9][10][11] , and perfect ideal squeezing is obtained at the SRPT critical point, as recently found both numerically and analytically 12,13 . In contrast to the standard squeezed state generation in non-equilibrium situations, quantum squeezing at the SRPT critical point is intrinsically stable and resilient against any noise even at finite temperatures 13 . As a result of this stable squeezing, such systems are intrinsically robust against decoherence, which is especially important for quantum sensing and continuous-variable quantum information technology.

A unique aspect of the SRPT is its manifestation as a physical phenomenon associated with the thermal-equilibrium state of a coupled light-matter system. This deviates from typical quantum-optics research that mainly deals with nonequilibrium excited-state dynamics. The occurrence of nonequilibrium SRPTs has been demonstrated in cold-atom systems driven by laser beams [14][15][16][17] . Although the temperature of cold atoms in a steady state is usually measured by the variance of their kinetic energy, the non-equilibrium SRPTs are inherently driven, dissipative, and transient phenomena. Effective temperatures defined with the driving power in the nonequilibrium SRPTs have been discussed theoretically 17 . However, the realization of SRPTs under pure thermal equilibrium is yet to be achieved. The existence of a SRPT analogue has been theoretically proposed for a superconducting circuit maintained under thermal equilibrium [18][19][20][21][22][23][24] , but no experimental observations of this effect have been reported.

The present work shows theoretically that the phase transition in erbium orthoferrite (ErFeO 3 ) with a critical temperature T c of ~4 K, known as the low-temperature phase transition (LTPT), is a magnonic SRPT, that is, an SRPT in which the Er 3þ spins cooperatively couple with the Fe 3þ magnonic field (spin-wave field) instead of with a photonic field as in the originally proposed SRPT. Specifically, we found that the LTPT occurs owing to Er 3þ -magnon coupling, even in the absence of direct Er 3þ -Er 3þ exchange interactions. In addition, we observed that the Er 3þ -magnon coupling enhances the T c value for LTPT compared to that obtained via direct Er 3þ -Er 3þ interactions. These results demonstrate the uniqueness of ErFeO 3 as a physical system in which SRPT can be experimentally realized under thermal equilibrium.

Principle of magnonic SRPT. The SRPT was first suggested in 1973 by Hepp and Lieb 1 , and has been extensively discussed based on the Dicke model 6 , conventionally expressed as

Here, â is the annihilation operator of a photon in a photonic mode with resonance frequency ω ph , Ŝx;y;z are spin N 2 operators representing an ensemble of two-level atoms with a transition frequency ω ex , and N is the number of atoms. The last term represents the coupling between the photonic mode and atomic ensemble with strength g. In the thermodynamic limit, i.e., in the limit of N !1, the SRPT arises when 4g 2 > ω ph ω ex , i.e., in the ultrastrong coupling regime g ≳ ω ph ; ω ex [3][4][5] . Below T c , the expectation values of the photon annihilation operator hâi and spin operator h Ŝz i become non-zero, indicating the spontaneous appearance of a static electromagnetic field and static polarization (or a persistent electric current) in thermal equilibrium.

The magnonic SRPT is a phase transition caused by ultrastrong coupling between a magnonic mode and other collective excitations in matter. The spontaneous appearance of magnons, also known as spin waves, reflects the spontaneous ordering of a spin ensemble mediating them in a certain direction. We present an explanation of the magnonic SRPT in the case of ErFeO 3 below.

Each unit cell in ErFeO 3 contains four Er 3þ ions and four Fe 3þ ions. The four Fe 3þ spins, each of which has an angular momentum of _S ¼ ð5=2Þ_, are oriented in different directions, even in the absence of an external direct current (DC) magnetic field 25 . However, it is known that the Fe 3þ spin resonances (magnon modes) may be described well by considering only two spins ŜA=B comprising two real Fe 3þ spins, which are usually treated as a single spin with S ¼ 5=2. In such a two-sublattice model of Fe 3þ , as depicted in Fig. 1a, the two spins ŜA=B are ordered antiferromagnetically along the c axis at T c < T ≲ 90 K, but are slightly canted toward the a axis and show weak magnetization (the Fe 3þ spins exhibit the so-called spinreorientation transition at 90 K ≲ T ≲ 100 K [26][27][28] ). In contrast, the Er 3þ spins are paramagnetic at T > T c , and they are directed along the a axis by the weak Fe 3þ magnetization. This phase is called the Γ 2 phase 29 .

At T < T c , as shown in Fig. 1b, when a two-sublattice model is used for the Er 3þ spins, they are ordered antiferromagnetically along the c axis, with a canting toward the a axis due to the Fe 3þ magnetization. Simultaneously, the Fe 3þ antiferromagnetism (AFM) vector S A À S B rotates gradually in the bc plane. The rotation angle measured from the c axis, φ, was estimated to be 49 at T ¼ 0 K 29 . This low-temperature phase is called the Γ 12 phase 29 .

The second-order phase transition between phases Γ 2 and Γ 12 at T c $ 4 K is called the LTPT 27,28 . There are at least two contributions to the LTPT, namely the Er 3þ -Er 3þ and Er 3þ -Fe 3þ exchange interactions 29,30 . Although the former is usually stronger than the latter and is largely responsible for LTPT, the latter is essential for explaining the rotation of the Fe 3þ AFM vector.

In the absence of the Er 3þ -Fe 3þ exchange interactions, as shown in Fig. 1a, the Fe 3þ spins are ordered antiferromagnetically along the c axis with a slight canting toward the a axis in the ground state of the Fe 3þ subsystem. Considering that the magnon excitation in this Fe 3þ subsystem corresponds to the photon excitation in the electromagnetic vacuum, the rotation of the Fe 3þ AFM vector (at T < T c as shown in Fig. 1b) indicates the spontaneous appearance of magnons, corresponding to the appearance of photons (a static electromagnetic field) in the ordinary SRPT, in thermal equilibrium. The ordering of the Er 3þ spins corresponds to the spontaneous appearance of an atomic field (polarization) in the SRPT.

Owing to the phenomenological similarities between the LTPT and SRPT, the microscopic models governing these phase transitions may be considered transferable. This constitutes the basic concept of magnonic SRPT proposed in this study.

ErFeO 3 spin model. First, we describe a general spin model for Er x Y 1Àx FeO 3 (0 ≤ x ≤ 1), which is consistent with our previous experimental study 31 . The replacement of Er 3þ ions with nonmagnetic Y 3þ ions simply reduces the density of the rare-earth (Er 3þ ) spins without changing the crystal structure or magnetic configuration of Fe 3þ spins in the Γ 2 phase 31,32 . Although the present work largely uses x ¼ 1 (ErFeO 3 ), the x-dependence is considered in "Spin resonance frequencies" in Supplementary Methods.

The Hamiltonian for the spins in Er x Y 1Àx FeO 3 consists of three parts, i.e.,

where H Fe , H Er , and H ErÀFe are the Hamiltonians of the Fe 3þ spins, Er 3þ spins, and Er 3þ -Fe 3þ interactions, respectively. As explained above, we employ the two-sublattice model for the Fe 3þ spins following Herrmann's model 33 and the methods of our prior works [31][32][33][34] . The Hamiltonian of the Fe 3þ spins is

Here, ŜA=B i is the operator of the Fe 3þ spin S ¼ 5=2 at the i-th site in the A/B sublattice, while ∑ n.n. represents a summation over all nearest-neighbor couplings. The number of nearest neighbors is

N 0 denotes the number of Fe 3þ spins in each sublattice and is equal to the unit cell count in ErFeO 3 . A total of 2N 0 spins represent the Fe 3þ subsystem. μ B is the Bohr magneton, and

is the g-factor tensor of the Fe 3þ spins. In the following, the g-factor of free electron spin is expressed as g. B DC is the external DC magnetic flux density. J Fe and D Fe y are the strengths of the isotropic and Dzyaloshinkii-Moriya-type exchange interaction strengths between the Fe 3þ spins, respectively. A x , A z , and A xz are the energies expressing the magnetic anisotropy of the Fe 3þ spins.

Although we expressed the Er 3þ subsystem using a single spin lattice for the paramagnetic Er 3þ spins (T > T c ) in our previous works 31,34 , in this study we employ a two-sublattice model for the Er 3þ spins to describe the Er 3þ -Er 3þ exchange interaction and LTPT. The Hamiltonian of the Er 3þ spins is

Here, RA=B i is the operator of rare-earth (Er 3þ or Y 3þ ) spin at the site i in the A/B sublattice. For Er x Y 1Àx FeO 3 , the rare-earth spins are represented randomly as s ¼ A; B; i.e.,

We describe each Er 3+

spin using a vector of Pauli operators σs i ðσ s i;x ; σs i;y ; σs i;z Þ t satisfying σs i;ξ σs i;ξ ¼ 1, ½σ s i;ξ ; σs 0 i 0 ;ξ ¼ 0 (ξ ¼ x; y; z), ½σ s i;x ; σs 0 i 0 ;y ¼ i2σ s i;z δ s;s 0 δ i;i 0 , ½σ s i;y ; σs 0 i 0 ;z ¼ i2σ s i;x δ s;s 0 δ i;i 0 , and ½σ s i;z ; σs 0 i 0 ;x ¼ i2σ s i;y δ s;s 0 δ i;i 0 , where δ i;j is the Kronecker delta. The Y 3þ ion is non-magnetic, and Rs i is replaced with 0. The first term in Eq. (6) represents the Zeeman effect, and the magnetic moment is expressed in terms of the anisotropic g factors g Er x;y;z for the Er 3þ spins as μs i À 1 2 μ B ðg Er x Rs i;x ; g Er y Rs i;y ; g Er z Rs

The factor 1=2 is added because ð1=2Þσ s i theoretically corresponds to a spin 1 2 operator. We define the g-factor tensor for the Er 3+ spins as

The second term in Eq. ( 6) represents the Er 3þ À Er 3þ exchange interaction with strength J Er . Because the Er 3þ ions are diluted in Er x Y 1Àx FeO 3 , the number of nearest-neighbor Er 3þ spins is effectively given by

We describe the Er 3þ -Fe 3þ exchange interactions as

In our model, the Er 3þ -Fe 3þ interaction is closed in each unit cell; that is, the Er 3þ and Fe 3þ spins in the same unit cell interact with each other but do not interact with the spins in other unit cells. J and D s;s 0 are the strengths of the isotropic and antisymmetric exchange interactions, respectively [31][32][33][34] . Considering the spin configuration at T < T c with no external DC magnetic field (see more details in "Reduction of number of parameters" in Supplementary Methods), we assume that D s;s 0 are expressed in terms of two values D x and D y as

As explained in detail in the section "Mean-field Calculation

Method" we assume that the y components RA=B i;y of the Er 3þ spins are not influenced by the Er 3þ -Fe 3þ interactions by implicitly considering a higher energy potential than that of the Er 3þ -Fe 3þ interaction strengths J and D s;s 0 along the b axis. This assumption helps to obtain an appropriate description of LTPT in accordance with our numerical calculations.

The actual values of the parameters appearing in our proposed spin model are provided in "Parameters" together with a description of how they were determined based on recent experimental results of terahertz magnetospectroscopy 31 and magnetization measurements 26 .

LTPT phase diagrams. Next, we show that our spin model certainly describes the thermal equilibrium (average) values of the Er 3þ spins σ A=B and Fe 3þ spins S A=B in the zero-wavenumber (infinite-wavelength) limit using the mean-field method. Details pertaining to the mean-field method are provided in the section "Mean-field Calculation Method." Because we simply considered a homogeneous external DC magnetic flux density B DC , σ A=B and S A=B were independent of the site index i. Figure 2a-c show the calculated phase diagrams as functions of T and B DC , applied along the a, b, and c axes, respectively. The difference j σ A z À σ B z j in the z components of the thermal equilibrium values of Er 3þ spins (AFM vector) is plotted in red. j σ A z À σ B z j is the order parameter for the LTPT in the presence of an external DC magnetic field in general, although the rotation angle of the Fe 3þ AFM vector can be utilized as an alternative order parameter if the external DC field is zero or is along the a axis. The bold solid curves represent phase boundaries.

These phase diagrams are consistent with those reported by Zhang et al. 26 . As shown in Fig. 2a, because ErFeO 3 possesses a weak magnetization along the a axis, the critical field depends on whether the field is parallel (in the same direction) or antiparallel (in the opposite direction) to the magnetization. The parameters used in the calculations are provided in "Parameters."

Figure 3 plots the thermal equilibrium values of the Er 3þ and Fe 3þ spins in the absence of an external DC magnetic field as functions of temperature. The LTPT, that is, the antiferromagnetic ordering of the Er 3þ spins along the c axis and the rotation of the Fe 3þ spins in the bc plane 29 , are reproduced well in our spin model. The rotation angle of the Fe 3þ AFM vector is φ ¼ 46 at T ¼ 0 K with our parameters. This value is approximately equal to the experimentally estimated value φ ¼ 49 29 .

Extended Dicke Hamiltonian. The mean-field method employed in this study, as illustrated in Figs. 2 and 3, is a standard means of analysing magnetic phase transitions. To investigate the analogy between LTPT and SRPT using the Dicke model, we derive an extended version of the Dicke model transformed from the spin model in Eq. ( 2). This derivation is given in detail in the section "Derivation of Extended Dicke Hamiltonian".

The extended Dicke Hamiltonian minimally including the terms relevant to the LTPT in an external DC magnetic field applied along the a axis, where the Γ 12 symmetry remains, is

-4 -3 -2 -1 0 1 2 3 0 1 2 3 4 0 0.5 1 1.5 2 0 2 4 6 0 1 2 3 4 0 0.5 1 1.5 0 1 2 3 4 External DC field B x (T) (a) (b) (c) External DC field B y (T) External DC field B z (T) Temperature T (K) Temperature T (K) Temperature T (K)

Fig. 2 Phase diagrams of spins in ErFeO 3 calculated using the mean-field method. An external DC magnetic field was applied along the a a. b b. c c axes. The difference j σ A z À σ B z j of the z components of the thermalequilibrium values of Er 3þ spins is mapped in red. The bold solid curves represent the phase boundaries. The external direct current (DC) magnetic field was varied from zero to positive or negative values at a fixed temperature. As ErFeO 3 shows weak magnetization along the a axis, the critical field depends on whether the field is parallel or antiparallel to the magnetization in Fig. 2a.

0 1 2 3 4 5 0 0.5 1 x y z 0 1 2 3 4 5 0 1 2 3 S x S y S z (a) Temperature T (K) Temperature T (K) Er 3+ spin Fe 3+ spin (b)

x 10

Fig. 3 Thermal equilibrium spin values. a Er 3þ spin. b Fe 3þ spin calculated using the mean-field method as functions of T in the case of zero external direct current (DC) magnetic field. As shown in Fig. 3a,

spontaneously appears below T c ¼ 4:0 K, i.e., the Er 3þ spins are antiferromagnetically ordered along the c axis. They show magnetization along the a axis as

due to the Er 3þ -Fe 3þ exchange interaction with the weak Fe 3þ magnetization along the a axis, whereas σ y ¼ σ

Here, âπ (â y π ) is the annihilation (creation) operator of an Fe 3þ magnon in the quasi-antiferromagnetic (qAFM) mode 33 . The eigenfrequency ω π ¼ 2π 0:896 THz can be evaluated using Eq. ( 63). The actual value was evaluated using the parameters shown in "Parameters." The Er 3þ resonance frequency is defined as follows.

The total number of 1 2 spins (Er 3þ spins) in the two sublattices is

Σx;y;z are spin N 2 operators representing the rare-earth spins (a detailed definition is given in Eq. ( 92)). The two Er 3þ -magnon coupling strengths in the last two terms of Eq. ( 16) are defined as follows.

Comparing Eq. ( 16) with Eq. ( 1) (the Dicke model), because âπ and Σx;y;z in Eq. ( 16) correspond to â and Ŝx;y;z in Eq. ( 1), respectively, we may observe that the g z term in Eq. ( 16) corresponds to the matter-photon coupling (transverse coupling), that is, the last term in Eq. ( 1). In addition, the g x term represents longitudinal coupling, and the term J Er describes the Er 3þ -Er 3þ exchange interactions in Eq. ( 16). The coupling strength g z ¼ 2π 0:116 THz shows the system fall into the ultrastrong regime because it is a considerable fraction of the Er 3þ resonance and qAFM magnon frequencies, E x ¼ h 0:023 THz and ω π ¼ 2π 0:896 THz. When the g z term causes an SRPT, h Σz i spontaneously acquires a non-zero value in thermal equilibrium, corresponding to the antiferromagnetic ordering of the Er 3þ spins along the c axis. As explained in "Derivation of Extended Dicke Hamiltonian," the spontaneous appearance of non-zero hiðâ y π À âπ Þi coupled with Σz in the g z term, corresponds to that of the Fe 3þ AFM vector in the b axis and causes its rotation in the bc plane. The Fe 3þ quasi-ferromagnetic (qFM) magnon mode can be neglected in describing the LTPT, because the AFM ordering of the Er 3þ spins and the spontaneous appearance of qAFM magnons are rather favoured and they prevent the appearance of qFM magnons, which feel additional energy cost under the ordering of the Er 3þ spins and the qAFM magnons.

As seen in Eqs. ( 20) and ( 21), the transverse coupling strength g z depends on D x , and the longitudinal coupling strength g x depends on J and D y . These expressions are reasonable from the perspective of the spin model in Eq. (11). The D x antisymmetric Er 3þ -Fe 3þ exchange interaction is essential for the LTPT because it couples σA=B z and ŜA=B y , which appear spontaneously at T < T c . In contrast, the J and D y exchange interactions are not directly related to the LTPT because these interaction terms do not couple σA=B z and ŜA=B y directly.

Evidence of magnonic SRPT. Using the semiclassical method described in "Semiclassical Calculation Method" with the extended Dicke Hamiltonian in Eq. ( 16), we calculated the thermal equilibrium values of the Er 3þ and Fe 3þ spins and magnon amplitudes as functions of temperature. Here and also in the calculation of the LTPT phase diagrams by the mean-field approach, we implicitly assumed that a thermal bath is connected to the extended Dicke Hamiltonian in Eq. (16) [and the spin model in Eq. ( 2)]. The thermal bath simply ensures that the system is in thermal equilibrium at a certain temperature in the present calculations, whereas it causes the energy loss and decoherence in non-equilibrium dynamics of Er 3þ spins and Fe 3þ magnons. Figure 4a-c show the thermal equilibrium values of the Er 3þ spins σ x;y;z ¼ h Σx;y;z i=ðN=2Þ, Fe 3þ spins S x;y;z , and Fe 3þ qAFM magnons a r;i as functions of temperature in the absence of an external DC magnetic field, i.e., with B DC ¼ 0. S x;y;z were calculated using Eqs. ( 113)-(115) with 4a, b, respectively, reproduce Fig. 3a, b calculated using the mean-field method with the original spin model, including T c , although S z differs. As depicted in Fig. 3b, a decrease in temperature causes a reduction in S z along with the spontaneous appearance of S y , whereas Fig. 4b reveals S z to remain nearly constant. This difference exists because

not hold in the extended Dicke Hamiltonian derived via magnon quantization (i.e., bosonisation of Fe 3þ spin modulations). The spins. c Fe 3þ magnon amplitudes as functions of T. These values were calculated using the semiclassical method with the extended Dicke Hamiltonian in the case of zero external direct current (DC) magnetic field. Figure 4a, b are nearly the same as Fig. 3a, b, respectively, except S z , which changes only slightly due to magnon bosonization. The Fe 3þ spins, S x;y;z , were calculated using Eqs. ( 113)-(115) with the thermal equilibrium value of the quasi-antiferromagnetic (qAFM) magnon annihilation operator

ultrastrong term g z , the last term in Eq. ( 16), causes the spontaneous appearance of σ z and a i , as shown in Fig. 4a, c, respectively, and the latter causes a non-zero S y through Eq. (114). The Fe 3þ AFM vector is rotated due to the spontaneous appearance of a non-zero S y when S

Thus, the LTPT, that is, the spontaneous ordering of Er 3þ spins (spontaneous appearance of σ z ) and the spontaneous rotation of Fe 3þ AFM vector (spontaneous appearance of a i and S y ), is caused by the Er 3þ -magnon coupling.

To compare the contributions of the Er 3þ -magnon couplings and Er 3þ -Er 3þ exchange interactions for the LTPT, Fig. 5 depicts the phase boundaries calculated using the full Hamiltonian (solid curves) as well as in the absence of Er 3þ -Fe 3þ exchange interactions (dashed-dotted curve;

and Er 3þ -Er 3þ exchange interactions (dashed curve; J Er ¼ 0). Figure 5a, b illustrate the results obtained using the mean-field and semiclassical methods with extended Dicke Hamiltonian, respectively. The solid curve in Fig. 5a is equal to that in Fig. 2a. The slight differences between Fig. 5a, b are discussed in "Aspects of phase boundaries" in Supplementary Methods.

The dashed curves (J Er ¼ 0) in Fig. 5 reveal that the phase transition occurs even in the absence of Er 3þ -Er 3þ exchange interactions and that T c equals approximately 1.2 K at B DC ¼ 0. Thus, Er 3þ -magnon coupling alone can cause the LTPT. In this sense, the LTPT can be interpreted as a magnonic SRPT because the Er 3þ -magnon coupling is sufficiently strong for the phase transition to occur.

On the other hand, in the absence of Er 3þ -magnon coupling, as denoted by the dashed-dotted curves, T c is approximately 2.6 K at B DC ¼ 0. This result appears to indicate that the contribution of the Er 3þ -Er 3þ exchange interactions is larger than that of the Er 3þ -magnon coupling. However, the actual T c is 4 K; that is, the Er 3þ -magnon coupling enhances the T c of the phase transition. In the same manner, the critical magnetic fields are also enhanced. These facts are similar to the suggestion of T c enhancement through photon-matter coupling by Mazza and Georges 35 ; however, in their case, phase transition does not occur solely by photon-matter coupling, and their model does not guarantee gauge invariance 36,37 .

Although the g z term causes the spontaneous appearance of both σ z and S y following the above-mentioned description of the SRPT, a non-zero σ z can also spontaneously appear due to the J Er term (Er 3þ -Er 3þ exchange interactions). Although Er 3þ -magnon coupling is inevitable for the spontaneous rotation of the Fe 3þ AFM vector (spontaneous appearance of S y ), we quantitatively evaluate the contributions of the Er 3þ -magnon coupling and Er 3þ -Er 3þ exchange interactions for the LTPT as follows.

The two contributions to the LTPT can be determined by analysing the condition for the SRPT in our extended Dicke Hamiltonian in Eq. ( 16) under the Holstein-Primakoff transformation [38][39][40] . The detailed calculations are discussed in the "Condition for SRPT in Extended Dicke Hamiltonian" section. The condition can finally be expressed as

For J Er ¼ g x ¼ 0, this expression is reduced to 4g z 2 > ω π ω Er for the SRPT in the Dicke model, Eq. ( 1).

The three terms on the left-hand side of Eq. ( 22) are evaluated as follows.

In the following, we refer to these quantities as coupling depths.

They are dimensionless measures of coupling strength and are determined based on the appearance of the SRPT. As seen in Eq. ( 22), the SRPT occurs when the sum of these coupling depths exceeds unity, i.e.,

The coupling depth D J Er of the J Er term is the largest, which is consistent with Fig. 5. The g x term (longitudinal coupling) has a negative contribution to the SRPT (D g x < 0). Among the three couplings, the contribution of the g z term is

23, and that of the total Er 3þ -magnon coupling is

These values are roughly equal to 1:3 K=ð1:3 K þ 3:4 KÞ ¼ 0:28, as estimated by Kadomtseva et al. 29 However, the longitudinal coupling (g x term) was not included in their model 41,42 , and the parameters were determined only by the phase boundary for B DC ==a.

Considering the analogy between an LTPT and SRPT, the coupling depth of the g z term satisfies D g z > 1 and

This result suggests that the transverse Er 3þ -magnon coupling is much stronger than the longitudinal coupling (giving a negative contribution) and is sufficiently strong to cause the SRPT alone. In this sense, we can conclude that the LTPT in ErFeO 3 is a magnonic SRPT obtained in the extended Dicke Hamiltonian with direct atom-atom interaction and longitudinal coupling (g x term).

As shown above, we quantitatively confirmed that the LTPT in ErFeO 3 is a magnonic version of the SRPT in thermal

-4 -3 -2 -1 0 1 2 0 1 2 3 4 5 -4 -3 -2 -1 0 1 2 0 1 2 3 4 5 External DC field B x (T) External DC field B x (T) Temperature T (K) Temperature T (K) Er 3+ -Fe 3+ & Er 3+ -Er 3+ Er 3+ -magnon & Er 3+ -Er 3+ Er 3+ -Er 3+ only Er 3+ -magnon only Er 3+ -Fe 3+ only Er 3+ -Er 3+ only 3 3 (a) (b)

Fig. 5 Phase boundaries of the low-temperature phase transition (LTPT) in ErFeO 3 . Boundaries calculated using the (a) mean-field method and (b) semiclassical method with the extended Dicke Hamiltonian. An external direct current (DC) magnetic field is applied along the a axis. The solid curves are the phase boundaries determined using the full Hamiltonian, and those in Figs. 5a and 2a are equivalent. The dashed-dotted curves are the phase boundaries in the absence of Er 3þ -magnon coupling (Er 3þ -Fe 3þ exchange interactions). The dashed curves are those obtained in the absence of Er 3þ -Er 3þ exchange interactions, i.e., the LTPT can be caused solely by the Er 3þ -magnon coupling and thus can be interpreted as a magnonic superradiant phase transition (SRPT).

equilibrium. This is the first confirmation of the SRPT since its proposal in 1973 1 .

Early reports on the SRPT suggested its no-go theorems [43][44][45][46] , implying that thermal-equilibrium SRPTs cannot be realized in systems described by the minimal-coupling Hamiltonian, that is, charged particles (without spins) interacting with electromagnetic fields. Because the classical treatment of the electromagnetic fields used in proofs of such no-go theorems can be justified only in limited situations 2,[45][46][47][48][49] , proposals of counter-examples against the no-go theorems and criticisms against the counter-examples have been repeated in SRPT research [35][36][37][50][51][52][53][54][55][56][57] .

One means of evading the no-go theorems involves introducing another degree of freedom, such as spin 44 . For example, it has been shown that the Rashba spin-orbit coupling can cause paramagnetic instability in an ultrastrongly coupled system between a cyclotron resonance and cavity photon field, implying an SRPT 37 . Further, it has been pointed out recently 37,58,59 that the coupling between matter and a spatially-varying multi-mode cavity fields plays a key role for circumventing the no-go theorem. Another method is to utilize various types of interactions and spin waves in magnetic materials, which cannot be described by the minimal-coupling Hamiltonian.

Ultrastrong photon-magnon coupling has been reported for an yttrium-iron-garnet sphere embedded in a cavity with a resonance frequency in the gigahertz region [60][61][62][63][64] , where the electromagnetic wave was confined by metallic or superconducting mirrors. Recently, g=ω $ 0:46 has been achieved to detect dark matter (galactic axions) 65 . Ultrastrong spin-magnon 31 and magnon-magnon 66,67 couplings have also been observed. However, evidence of an SRPT has not been reported even with those magnonic ultrastrong couplings, although various phase transitions exist in magnetic systems, and it is conceivable that some of the known phase transitions can be understood as the SRPT or an analogue.

The LTPT in ErFeO 3 has been discussed in relation to the cooperative Jahn-Teller transition 29,41,42 , which is analogous to the SRPT 68,69 . Vitebskii and Yablonskii proposed a theoretical model for describing the LTPT in 1978 30 . Further, Kadomtseva et al. theoretically investigated the ratio between the Er 3þ -Er 3þ and Er 3þ -Fe 3þ interaction strengths in 1980 29 . They also mentioned the analogy between the LTPT and cooperative Jahn-Teller transition 41,42 . Loos and Larson discussed the analogy between the cooperative Jahn-Teller transition and SRPT in 1984 and 2008, respectively 68,69 . However, the analogy between the LTPT and SRPT has not been directly drawn either theoretically or experimentally because the analogies between the LTPT and cooperative Jahn-Teller transition and between the latter and the SRPT have been independently discussed 29,68,69 , and no experimental evidence has been shown. The spin-Peierls transitions 70,71 and the spin-reorientation transition in rare-earth iron garnets 72 have also been discussed as analogous phenomena to the SRPT. However, no experimental evidence has been demonstrated. Structural transitions in ferroelectric materials may also be seen as a SRPT analogue at a first glance 73 . However, when we map such ferroelectric systems to the Dicke model, we find that the resonance frequency of the electric polarization becomes an imaginary value, which indicates that such ferroelectric phase transitions are caused by the instability of the electric polarization subsystem rather than by the coupling between the polarization and phonon subsystems. Hence, no other SRPT analogue by matter-matter coupling has been confirmed quantitatively.

In 2018, the ffiffiffiffi N p -dependence (N is the Er 3þ density) of the anticrossing frequency, or vacuum Rabi splitting (2g), between paramagnetic Er 3þ spins and a Fe 3þ magnon mode was experimentally confirmed at T > T c 31 . This ffiffiffiffi N p -dependence, the Dicke cooperativity, can be taken as evidence that the coupling between the Er 3þ spin ensemble and Fe 3þ magnon mode is cooperative, being well described by the Dicke model or its extension.

This study provides the quantitative evidence of magnonic SRPT manifestation. Meanwhile, the existence of photonic SRPT, which was originally proposed in 1973, has yet to be confirmed. Moreover, the possibility of its theoretical existence in materials with spin degree of freedom is still under debate 36,37,58,59 . Since the development of the Dicke model, this study is the first to elucidate the occurrence of magnonic thermal-equilibrium SRPT in an actual material, namely ErFeO 3 .

In this study, using an ErFeO 3 spin model reproducing both the phase diagrams obtained via magnetization measurements 26 and terahertz magnetospectroscopy results 31 , we derived an extended Dicke Hamiltonian that accounts for Er 3þ -Er 3þ exchange interactions as well as the cooperative coupling between the Er 3þ spins and Fe 3þ magnon modes. We found that the LTPT in ErFeO 3 can be caused solely by Er 3þ -magnon coupling (in the absence of Er 3þ -Er 3þ exchange interactions). From the analytical correspondence between the spin and Dicke models and the quantitative verification that the Er 3þ -magnon coupling solely causes the LTPT, we concluded that the LTPT in ErFeO 3 is a magnonic SRPT in the extended Dicke model. This is the first confirmation of the SRPT in thermal equilibrium since its proposal in 1973 1 . These results are expected to be the first step in finding the (originally proposed) photonic SRPT in magnetic or other materials explicitly including the spin degree of freedom.

The thermal SRPT in ErFeO 3 would exhibit rich physics beyond the quantum or zero-temperature SRPT that has been demonstrated in laser-driven cold atoms [14][15][16][17] . It is known that the thermal and quantum fluctuations of photons and atoms exhibit characteristic behaviours around the SRPT 74,75 . Recent studies have reported the occurrence of strong, two-mode quantum squeezing at the SRPT critical point 12,13 . In future endeavours, including ongoing terahertz magnetospectroscopy experiments on Er x Y 1Àx FeO 3 concerning LTPT 76 and subsequent quantumfluctuation measurements 77,78 of magnons and Er 3+ spins, we intend to investigate the occurrence of such quantum-squeezing phenomena during thermal SRPT.

The generation of squeezed states of light has attracted considerable research interest over several decades because they facilitate precision measurements to be performed beyond the limitations encountered owing to the manifestation of quantum vacuum fluctuations and evolution of continuous-variable quantum computing. However, most existing squeezing-generation protocols require a system to be driven to realize transient squeezed states. This limits the realizable degree of squeezing during experiments owing to unpredictable noise. In contrast, quantum squeezing at the SRPT critical point can be stably realized under thermal equilibrium, because an ultrastrong coupled system remains at its most stable in the squeezed state. As a result, such systems are resilient to unpredictable noise. This fundamental stability and resilience are expected to facilitate the realization of novel applications exploring quantum sensing and decoherence-robust continuous-variable quantum computing via the occurrence of quantum squeezing at the SRPT critical point.

Parameters. Following our previous study 31 , we used the following values for the Fe 3þ subsystem in our numerical calculations, except A x , which was determined to fit the spin resonance frequencies to the corresponding terahertz absorption spectrum in our experiments 31 (see "Spin resonance frequencies" in Supplementary Methods for details).

The anisotropic g-factors for Er 3þ spins were assumed to be

These values were utilized to fit the Er 3þ spin resonance frequencies depicted in Supplementary Figs. 123to their corresponding absorption peak positions observed during experiments 31 (refer to "Spin resonance frequencies" in Supplementary Methods). They were multiplied by 2 compared to those estimated in our previous study 31 to compensate for the use of the additional factor of 1=2 in Eq. ( 8). The anisotropic g-factors for Fe 3þ spins were assumed to be

Here, g Fe z was determined to reproduce the critical magnetic flux density B DC z $ 20 T 26 of the transition between the Γ 2 phase and the Γ 4 phase, in which the Fe 3þ spins are ordered antiferromagnetically along the a axis with slight canting toward the c axis, in the case of B DC ==c. On the other hand, g Fe x and g Fe y were simply set to the values in the case of free electron spin because the results in the present study are insensitive to these values.

Concerning the Er 3þ -Er 3þ and Er 3þ -Fe 3þ exchange interactions, we used the following values.

These values were utilized to fit Fig. 2 roughly to the phase diagrams reported by Zhang et al. 26 . The precise values of J Er , J, and D y were mainly determined to fit our calculated spin resonance frequencies B DC ==c to the corresponding terahertz absorption spectrum in our experiments 31 , which are both shown in Supplementary Fig. 3a (see "Spin resonance frequencies" in Supplementary Methods). On the other hand, D x was determined to reproduce T c ¼ 4:0 K. Although the ratio between the Er 3þ ÀEr 3þ and Er 3þ -Fe 3þ interaction strengths was theoretically investigated by the phase boundary for B DC ==a 29 , the phase diagrams (T c and critical DC fields) themselves were not sufficient to determine all of our parameters, although we do not intend to claim the impossibility of such determination in the present study. The phase diagrams gave only some ranges of the parameters. Because the LTPT is caused by not only the Er 3þ -Er 3þ exchange interaction, but also the Er 3þ -Fe 3þ interactions (Er 3þ -magnon couplings), there are at least four parameters (J Er , J, D x , and D y ) even if the number of parameters is reduced according to the analysis in "Reduction of number of parameters" in Supplementary Methods. The anisotropic g-factors g Er x , g Er y , and g Er z of the Er 3þ spins are free parameters, and could easily change the critical DC fields. T c and three critical DC fields obtained from the magnetization measurements 26 were not sufficient to determine the above parameters.

In determining all of these quantities, the spin resonance frequencies were informative. In particular, as discussed in "Spin resonance frequencies" in Supplementary Methods using the extended Dicke Hamiltonian, the Er 3þ -Er 3þ exchange interaction strength J Er clearly appears as the frequency splitting between the Er 3þ in-phase and out-of-phase resonances. The out-of-phase mode cannot be excited by the terahertz wave unless it couples with the Fe 3þ magnon modes. In that sense, the anti-crossing between the Er 3þ in-phase resonances, out-of-phase resonances, and Fe 3þ qFM magnon mode B DC z $ 4 T in Supplementary Fig. 3 provides the most important information for determining J Er and the other parameters (see "Spin resonance frequencies" in Supplementary Methods).

Mean-field calculation method. Because we simply considered a homogeneous B DC in this study, the expectation values of the Er 3þ spins σ A=B hσ A=B i i and Fe 3þ spins S A=B h ŜA=B i i were independent of the site index i. The brackets represent the theoretical expectation values of the operators at a finite temperature in the Heisenberg representation. The brackets also correspond to the ensemble average of the spins in each sublattice. Their equations of motion can be obtained from the Heisenberg equations derived from the Hamiltonian in Eq. ( 2), as follows (s ¼ A; B).

Here, B A=B Er and B A=B Fe are the mean fields for the Er 3þ and Fe 3þ spins, respectively, and they can be expressed as

In Eqs. ( 43) and ( 44), the first, second, and third terms represent the Zeeman effect, Er 3þ ÀEr 3þ exchange interaction, and Er 3þ ÀFe 3þ exchange interaction, respectively. In Eqs. ( 45) and ( 46), the first, second, and third terms represent the Zeeman effect, Er 3þ ÀFe 3þ exchange interaction, and Fe 3þ ÀFe 3þ exchange interaction, respectively. The dilution of the Er 3þ spins is reflected by the factors z Er ¼ 6x and x. z Er denotes the number of neighbours of Er 3þ , and its value effectively decreases by a factor of x. Because ð1=2Þσ A=B corresponds to the spin 1 2 operator, a factor of 2 appears overall in Eqs. (43) and (44). As explained at the end of the section "ErFeO 3 ," the y component of the third term in Eqs. ( 43) and ( 44) is set to zero by means of implicitly considering a high-energy potential.

The free energy of the system is minimized when the thermal equilibrium values (time averages) of spins σ A=B and S A=B are parallel to their mean fields B s Er B s Er ðf σ A=B g; f S A=B gÞ and B s Fe B s Fe ðf σ A=B g; f S A=B gÞ as follows.

Here, the unit vectors of the mean fields are defined as

The thermal equilibrium values σ A=B and S A=B can be determined as follows.

For given mean fields B s

s Er , the effective Hamiltonians of each Er 3þ and Fe 3þ can be defined as

Subsequently, the partition functions can be expressed as

Er Tr e À Ĥs Er =ðk B TÞ

where the following are defined.

Because σA=B is not a standard spin operator with an angular momentum of _ or _=2 but is a vector of the Pauli operators, the summation is performed for m ¼ ±1.

The free energies are given as Àk B TlnZ A=B Er and Àk B TlnZ A=B Fe , and the thermal equilibrium values of the spins are

where B S ðzÞ is the Brillouin function, defined as

By consistently solving Eqs. ( 43)-( 48), ( 57) and ( 58), σ A=B and S A=B can be determined at finite temperatures.

Derivation of extended Dicke Hamiltonian. Here, we describe the transformation of our spin model, Eq. ( 2), into an extended version of the Dicke Hamiltonian, Eq. ( 1). We first rewrite the Fe 3+ subsystem H Fe in terms of the annihilation and creation operators of a magnon in "Fe 3þ subsystem." The Er 3þ subsystem H Er is rewritten using large spin operators in "Er 3þ subsystems." The Er 3þ -Fe 3þ exchange interactions, H ErÀFe , are transformed into five Er 3þ -magnon couplings as per "Er 3þ ÀFe 3þ interactions." Finally, the extended Dicke Hamiltonian is discussed in the section "Total Hamiltonian." Fe 3þ subsystem. We assume that the most stable values of the Fe 3þ spins at zero temperature, S A=B , are unchanged even when B DC (≲10 T) is applied, as we also assumed in our previous studies 31,34 . Under this assumption, as depicted in Fig. 1a, the most stable state (i.e., the ground state) of the Fe 3þ subsystem ĤFe , Eq. (3), can be expressed as

Here, the canting angle β 0 can be expressed as (see "Magnon quantization" in Supplementary Methods and Supplementary Fig. 4 or refs. 31,33,34 )

The magnon is the quantum of spin fluctuations (spin waves) from this stable state. As shown in "Magnon quantization" in Supplementary Methods as well as in refs. 31,34 in the long-wavelength limit, the Fe 3þ Hamiltonian ĤFe , Eq. ( 3), can be rewritten in terms of the annihilation (creation) operators âK (â y K ) of Fe 3þ magnons as

Here, K ¼ 0 and π correspond to the qFM and qAFM magnon modes, respectively 33 . The eigenfrequencies can be expressed as follows.

Here, we define

The operators of the spin fluctuations δ ŜA=B

can be expressed as

where the following are defined.

For the subsequent discussion, we define the sum and difference of the spins as

Their equilibrium (most stable) values are

and their fluctuations are given by the sum and difference of Eqs. ( 68) and ( 69) as follows.

Er 3þ subsystems. We define the following new operators.

For an Er 3þ ion, ð1=2Þ RA=B i is a spin 1 2 operator and ΣA=B is a spin N 4 operator representing the rare earth spins in the A/B sublattice. We also define the sum and difference of the two sublattice spins as

In the long-wavelength limit, all spins in each sublattice have the same values in both static and dynamic situations. Subsequently, the Er 3þ Hamiltonian in Eq. ( 6) can be rewritten as

Er 3þ -Fe 3þ interactions. In the same manner as in our prior works 31,34 , the Hamiltonian of the Er 3þ -Fe 3þ exchange interactions can be rewritten using Eq. ( 11), as

In each set of parentheses, the first term represents the influence of the static components (equilibrium values) S A=B 0 of Fe 3þ spins to Er 3þ spins Σ ± , and the second term represents the coupling between the Fe 3þ fluctuation δ Ŝ ± and Er 3þ spins Σ ± . We divide these terms into two Hamiltonians as The first term gives part of the Er 3þ spin resonance frequency and can be expressed as follows.

Here, we used Eqs. ( 73), ( 74) and (18). We neglected ðÀ4SD x cos β 0 Þ ΣÀ y under the assumption explained at the end of the section "ErFeO 3 ." The second term in Eq. ( 82) can be rewritten in terms of the Fe 3þ fluctuations as

Total Hamiltonian. In terms of the annihilation and creation operators of magnons, the total Hamiltonian can be expressed as

The five coupling strengths are defined as

The actual values were evaluated using the parameters shown in "Parameters."

Compared with the expressions in our previous studies 31,34 , the coupling strengths in Eqs. ( 86)-(87) include additional factors ffiffi ffi 2 p and ffiffi ffi S p . First, ffiffi ffi 2 p

originates from the number of Er 3þ sublattices in the present study, whereas a single Er 3þ lattice was considered in our previous studies 31,34 . The second factor, ffiffi ffi S p , is a result of the difference in the method of normalizing the Fe 3þ spins between the present and previous studies 31,34 .

Whereas the Er 3þ spin ensemble is described by six operators Σþ x;y;z and ΣÀ x;y;z in the extended Dicke Hamiltonian in Eq. (85), only Σþ x and ΣÀ z are relevant to the LTPT shown in Fig. 1.

x corresponds to the paramagnetic alignment by the Fe 3þ magnetization along the a axis, and ΣÀ z corresponds to the antiferromagnetic ordering along the c axis. Subsequently, to analyse the thermal equilibrium values of the spins, it is sufficient to consider only the following two terms in the Er 3þ -Er 3þ exchange interactions.

In contrast, while the Fe 3þ spins are described by the qFM and qAFM magnon modes in Eq. ( 85), only the qAFM mode is relevant to the LTPT. As shown in Fig. 1, δ ŜÀ y and δ ŜÀ z are required to describe the rotation of the Fe 3þ AFM vector in the bc plane, and δ Ŝþ x is required for possible modulation of canting along the a axis. As seen in Eqs. ( 75) and ( 76), δ Ŝþ x , δ ŜÀ y , and δ ŜÀ z are related to the qAFM magnon mode (K ¼ π), and the qFM mode (K ¼ 0) plays no role in the LTPT.

Consequently, among the terms in the total Hamiltonian given by Eq. (85), it is only necessary to consider the terms shown in Eq. ( 16) to describe the LTPT (the other terms are required to fully reproduce the terahertz spectra discussed in "Spin resonance frequencies" in Supplementary Methods). Note that, in Eq. ( 16), we rewrote the large spin operators representing the Er 3þ spin ensemble as

where we re-indexed the Pauli operators representing the Er 3þ spins in the two sublattices as σA i;x ! σ2iÀ1;x σA i;y ! σ2iÀ1;y σA i;z ! σ2iÀ1;z 8 > > < > > : ; σB i;x ! σ2i;x σB i;y ! Àσ 2i;y σB i;z ! Àσ 2i;z 8 > > < > > : :

Further, in Eq. ( 16), it was assumed that the external DC magnetic field is applied along the a axis to maintain Γ 12 symmetry, where either j σ A z À σ B z j or the rotation angle φ of the Fe 3þ AFM vector from the c axis can be the order parameter for the LTPT. Among the five Er 3þ -magnon couplings in Eq. ( 85), only the g x and g z terms are required to consider the coupling between Σx;z and the qAFM magnons. Although the g y 0 term also couples Σy and qAFM magnons, its coupling strength is negligible compared with g x;z , as shown in Eq. ( 88), which is consistent with the experimentally observed antiferromagnetic ordering of the Er 3þ spins along the c axis (h ΣÀ y i ¼ 0). As demonstrated in Figs. 4 and 5, the LTPT can be quantitatively reproduced as the SRPT in the extended Dicke Hamiltonian, Eq. ( 85), which was derived from the spin model of ErFeO 3 . The essential terms were extracted as shown in Eq. ( 16). The g z term (antisymmetric Er 3þ -Fe 3þ exchange interaction with D x ) corresponds to the matter-photon coupling and causes the antiferromagnetic ordering of Er 3þ spins along the c axis and the b component of the Fe 3þ spins through the spontaneous appearance of qAFM magnons.

Semiclassical calculation method. Wang and Hioe demonstrated a simple calculation method for the SRPT in 1973 2 , and Hepp and Lieb confirmed its validity for the Dicke model 47 . In the Dicke model, the partition function at temperature T Z Dicke ðTÞ Tr½e À ĤDicke =ðk B TÞ ; ð94Þ in the thermodynamic limit N ! 1 can be approximately evaluated by replacing the trace over the photonic variables with an integral over coherent states j ffiffiffiffi N p a ( a 2 C; giving âj ffiffiffiffi N p a ¼ ffiffiffiffi N p aj ffiffiffiffi N p a ) as Z Dicke ðTÞ Z d 2 a π=N Tr½e À Ĥeff Dicke ð aÞ=ðk B TÞ ð95Þ ¼ Z d 2 a π=N e À S Dicke ð a;TÞ=ðk B TÞ ;

where an effective Hamiltonian is defined as given below,

as well as an action

and an effective Hamiltonian per atom

The normalized expectation value a ¼ hâi= ffiffiffiffi N p of the annihilation operator of a photon at temperature T can be determined to minimize the action, that is, ∂ S=∂Re½ a ¼ 0 and ∂ S=∂Im½ a ¼ 0. a acquires a non-zero value below T c when 4g 2 > ω ph ω ex is satisfied ( ffiffiffiffi N p a gives a finite electric (displacement) field or vector potential even in the thermodynamic limit N ! 1 if the atomic density is fixed). The above approximation is justified if the free energy F Dicke ðTÞ Àðk B T=NÞln Z Dicke ðTÞ per atom satisfies _ω ph =N ( j F Dicke ðTÞj in the thermodynamic limit 45,46,48,49 . Following the above treatment, we calculated the expectation values of the Er 3þ spin and Fe 3þ qAFM magnon operators in the extended Dicke Hamiltonian, Eq. ( 16), at a finite temperature. In the thermodynamic limit N ! 1 the partition function ZðTÞ Tr½e À Ĥ=ðk B TÞ can be approximately evaluated by replacing the trace over the magnonic variables with an integral over c-numbers a r ; a i 2 R,

, as ZðTÞ Z d a r d a i π=N Tr e À Ĥeff ð a r ; a i Þ=ðk B TÞ h i ð101Þ ¼ Z d a r d a i π=N e À Sð a;TÞ=ðk B TÞ ;

where we define an effective Hamiltonian Ĥeff ð a r ;

by introducing the Er 3þ components h Σx;z i of the mean fields of the Er 3þ ensemble. The action in Eq. ( 102) is defined as follows.

Sð a r ; a i ; TÞ Àk B TlnTr e À Ĥeff ð a r ; a i Þ=ðk B TÞ

Here, we define an effective Hamiltonian per Er 3þ spin as

The site index i is omitted here because all the spins are identical. The action S is minimized at ∂ S=∂ a r ¼ 0 and ∂ S=∂ a i ¼ 0, yielding

where the expectation values of the Pauli operators are defined for a given a r and a i , as σ ξ hσ ξ i Tr½σ ξ e À Ĥa ð a r ; a i Þ=ðk B TÞ Tr½e À Ĥa ð a r ; a i Þ=ðk B TÞ : ð109Þ

From Eqs. (107) and (108), the expectation values of the large spin operators can be expressed as

Substituting these into Eq. ( 106) gives

By simultaneously solving Eqs. ( 107)-( 109) and (112) for a given temperature, T, we obtain the thermal equilibrium values of the Er 3þ spins σ x;z and Fe 3þ qAFM magnons a r;i . From Eqs. ( 60) and ( 68)-( 71), the thermal equilibrium values of the Fe 3þ spins can be obtained from those of the qAFM magnons a r;i as

Condition for SRPT in extended Dicke Hamiltonian. To quantitatively evaluate the contributions of the Er 3þ -magnon couplings and Er 3þ -Er 3þ exchange interactions to the LTPT, the condition for the SRPT in our extended Dicke Hamiltonian, Eq. ( 16), can be derived using the Holstein-Primakoff transformation [38][39][40] . Σx;y;z can be rewritten using the bosonic annihilation (creation) operator b ( by ) as Σx ! by b À N 2 ; ð116Þ Σy ! by ðN À by bÞ 1=2 þ ðN À by bÞ 1=2 b 2 ; ð117Þ Σz ! by ðN À by bÞ 1=2 À ðN À by bÞ 1=2 b i2 :

Further, all the operators are replaced by c-numbers a r ; a i ; b 2 R as

Subsequently, the Hamiltonian in Eq. ( 16) becomes

The ground state of the system should satisfy

Solving the first two equations, the Fe 3þ qAFM magnon amplitudes can be expressed as

By substituting these expressions into Eq. ( 124), the following equation for the Er 3þ amplitude can be obtained.

For a real non-zero value of b to exist, the parameters must satisfy Eq. ( 22).