Terahertz Strong-Field Physics without a Strong External Terahertz Field

TL;DR

This paper explores how matter can exhibit strong-field phenomena through quantum vacuum interactions in cavity systems, especially in the ultrastrong coupling regime, without the need for external terahertz fields.

Contribution

It discusses two condensed matter systems demonstrating cooperative ultrastrong interactions in the terahertz range, expanding understanding of vacuum-field coupled systems.

Findings

Observation of vacuum Rabi splitting in systems without external fields

Identification of systems reaching the ultrastrong coupling regime

Potential for quantum phase transitions in matter-vacuum interactions

Abstract

Traditionally, strong-field physics explores phenomena in matter (atoms, molecules, and solids) driven by an extremely strong laser field nonperturbatively. However, even in the complete absence of an external electromagnetic field, strong-field phenomena can arise when matter strongly couples with the zero-point field of the quantum vacuum state, i.e., fluctuating electromagnetic waves whose expectation value is zero. Some of the most striking examples of this occur in a cavity setting, in which an ensemble of two-level atoms resonantly interacts with a single photonic mode of vacuum fields, producing vacuum Rabi splitting. In particular, the nature of the matter-vacuum-field coupled system fundamentally changes when the coupling rate (equal to one half of the vacuum Rabi splitting) becomes comparable to, or larger than, the resonance frequency. In this so-called ultrastrong coupling regime, a non-negligible number of photons exist in the ground state of the coupled system. Furthermore, the coupling rate can be cooperatively enhanced (via so-called Dicke cooperativity) when the matter is comprised of a large number of identical two-level particles, and a quantum phase transition is predicted to occur as the coupling rate reaches a critical value. Low-energy electronic or magnetic transitions in many-body condensed matter systems with large dipole moments are ideally suited for searching for these predicted phenomena. Here, we discuss two condensed matter systems that have shown cooperative ultrastrong interactions in the terahertz frequency range: a Landau-quantized two-dimensional electron gas interacting with high-quality-factor cavity photons, and an Er$^{3+}$ spin ensemble interacting with Fe$^{3+}$ magnons in ErFeO$_3$.

Figures

Equations

$${\hat{\mathcal{H}}}_{\text{Tavis–Cummings}}$$

$${\hat{\mathcal{H}}}_{\text{min}}=\int \mathrm{d}\mathbf{r}\left[\frac{{\epsilon }_0{\hat{\mathbf{E}}}_{\perp }(\mathbf{r}{)}^2}{2}+\frac{\hat{\mathbf{B}}(\mathbf{r}{)}^2}{2{\mu }_0}\right]+\sum _{i=1}^N\left[\frac{[{\hat{\mathbf{p}}}_i+e\hat{\mathbf{A}}({\hat{\mathbf{r}}}_i){]}^2}{2m}+V(\{{\hat{\mathbf{r}}}_i\})\right].$$

$${\hat{\mathcal{H}}}_{\text{Dicke}}$$

$${\hat{\mathcal{H}}}_{\text{Hopfield}}$$

$${\hat{\mathcal{H}}}_{\text{Landau}}$$

$${\hat{\mathcal{H}}}_{\mathrm{ErFe}{\mathrm{O}}_3}$$

$${\hat{\mathcal{H}}}_{\text{Landau}}={\hat{\mathcal{H}}}_{\text{EM}}+{\hat{\mathcal{H}}}_{\text{2DEG}}.$$

$${\hat{\mathcal{H}}}_{\text{EM}}$$

$$=\sum _{\xi =x,y}\int \mathrm{d}z\left[\frac{{\hat{\Pi }}_{\xi }(z{)}^2}{2{\epsilon }_0{\epsilon }_{\text{cav}}(z)}+\frac{1}{2{\mu }_0}{\left(\frac{\partial }{\partial z}{\hat{A}}_{\xi }(z)\right)}^2\right].$$

$${\hat{\mathcal{H}}}_{\text{2DEG}}$$

$$\frac{{\partial }^2}{\partial z^2}{\hat{A}}_{\xi }(z,t)-\frac{{\epsilon }_{\text{cav}}(z)}{c^2}\frac{{\partial }^2}{\partial t^2}{\hat{A}}_{\xi }(z,t)=0.$$

$$\frac{1}{\sqrt{{\epsilon }_{\text{cav}}(z)}}\frac{{\partial }^2}{\partial z^2}\frac{f_k(z)}{\sqrt{{\epsilon }_{\text{cav}}(z)}}+\frac{{\omega }_k{}^2}{c^2}{f}_k(z)=0.$$

$${\hat{A}}_{\xi }(z)$$

$${\hat{\Pi }}_{\xi }(z)$$

$${\hat{\mathcal{H}}}_{\text{EM}}=\sum _k\sum _{\xi =x,y}\hslash {\omega }_k\left({\hat{a}}_{k,\xi }^{†}{\hat{a}}_{k,\xi }+\frac{1}{2}\right).$$

$${\hat{\mathcal{H}}}_{\text{2DEG}}$$

$${\hat{A}}_{\pm }(z)$$

$$=\sum _k\sqrt{\frac{\hslash }{2{\epsilon }_0{\epsilon }_{\text{cav}}(z){\omega }_{k}S}}{f}_k(z)\left({\hat{a}}_{k,\mp }^{†}+{\hat{a}}_{k,\pm }\right),$$

$${\hat{a}}_{k,\pm }\equiv \frac{{\hat{a}}_{k,x}\mp \mathrm{i}{\hat{a}}_{k,y}}{\sqrt{2}}.$$

$$\hat{b}\equiv \frac{1}{\sqrt{N}}\sum _{i=1}^N{\hat{c}}_i,$$

$$\mathrm{i}\sqrt{\frac{\hslash {\omega }_{\text{cyc}}{e}^2}{m}}\sum _{i=1}^N\left[{\hat{c}}_i^{†}{\hat{A}}_{+}(z_{\text{2DEG}})-{\hat{A}}_{-}(z_{\text{2DEG}}){\hat{c}}_j\right]$$

$$=\sum _k\mathrm{i}\hslash g_k\left[{\hat{b}}^{†}({\hat{a}}_{k,+}+{\hat{a}}_{k,-}^{†})-({\hat{a}}_{k,-}+{\hat{a}}_{k,+}^{†})\hat{b}\right],$$

$$g_k\equiv \sqrt{\frac{e^2{\omega }_{\text{cyc}}{n}_{\text{2DEG}}}{2{\epsilon }_0{\epsilon }_{\text{cav}}(z_{\text{2DEG}})m{\omega }_k}}{f}_k(z_{\text{2DEG}}).$$

$$\frac{N{e}^2}{2m}{\hat{\mathbf{A}}}^2=\sum _{k,k^{\prime }}\frac{\hslash g_k{g}_{k^{\prime }}}{{\omega }_{\text{cyc}}}({\hat{a}}_{k,-}+{\hat{a}}_{k,+}^{†})({\hat{a}}_{k^{\prime },+}+{\hat{a}}_{k^{\prime },-}^{†}).$$

$$\frac{{\omega }_a{}^2}{{\omega }_{\pm }{}^2}$$

$$=1-\frac{2{\tilde{g}}^2}{{\omega }_{\pm }({\omega }_{\pm }\mp {\omega }_{\text{cyc}})}.$$

$$\frac{{\omega }_a{}^2}{{\omega }^2}=1+\frac{4\tilde{g}{}^2}{{\omega }_b{}^2-{\omega }^2},$$

$$\hat{\mathcal{H}}={\hat{\mathcal{H}}}_{\text{Fe}}+{\hat{\mathcal{H}}}_{\text{Er}}+{\hat{\mathcal{H}}}_{\text{Fe–Er}}.$$

$${\hat{\mathcal{H}}}_{\text{Fe}}$$

$$-\sum _{i=1}^N\left(A_x{\hat{S}}_{i,x}^B{}^2+A_z{\hat{S}}_{i,z}^B{}^2-A_{xz}{\hat{S}}_{i,x}^B{\hat{S}}_{i,z}^B\right).$$

$${\hat{\mathcal{H}}}_{\text{Er}}=-\sum _{i=1}^N{\hat{\mathbf{\mu }}}_i\cdot {\mathbf{B}}_{\text{stat}}.$$

$${\hat{\mathcal{H}}}_{\text{Fe–Er}}=\sum _{i=1}^N\left[J_A{\hat{\mathbf{R}}}_i\cdot {\hat{\mathbf{S}}}_i^A+J_B{\hat{\mathbf{R}}}_i\cdot {\hat{\mathbf{S}}}_i^B+{\mathbf{D}}_A\cdot ({\hat{\mathbf{R}}}_i\times {\hat{\mathbf{S}}}_i^A)+{\mathbf{D}}_B\cdot ({\hat{\mathbf{R}}}_i\times {\hat{\mathbf{S}}}_i^B)\right].$$

$$\delta {\hat{\mathbf{S}}}_i^A$$

Full text

I Introduction

Matter in the presence of a strong electromagnetic (EM) field exhibits intriguing phenomena that cannot be understood by treating the field as a small perturbation Bucksbaum (1990). For example, the AC Stark effect, which occurs in the presence of a resonant laser field, represents coupled light-matter states; a measure of the coupling strength is the Rabi energy, $\hbar\Omega_{\text{Rabi}}=d_{12}E$, where $d_{12}$ is the transition dipole moment and $E$ is the field strength. In the nonresonant case, laser-driven matter exhibits other nonperturbative phenomena (e.g., above-threshold ionization and high-order harmonic generation), in which the ponderomotive potential $U_{\text{p}}=e^{2}E^{2}/(4m\omega^{2})$, i.e., the cycle-averaged quiver energy of an electron in laser light with frequency $\omega$ Bucksbaum (1990); Chin et al. (2000); Srivastava et al. (2004); Wegener (2005), plays a critical role. In both the resonant and nonresonant cases, different regimes of strong-field phenomena arise, depending on the normalized parameters $\hbar\Omega_{\text{Rabi}}/\hbar\omega$ and $U_{\text{p}}/\hbar\omega$, respectively. Since these parameters scale as $1/\omega$ and $1/\omega^{3}$, respectively, for a given field strength $E$, it is inherently advantageous to use smaller laser frequencies to explore strong-field physics. In particular, in condensed matter systems in the resonant case, there exist transitions with enormous $d_{12}$ in the terahertz (THz) frequency range, which makes it possible to explore uncharted regimes of strong-field physics even without a strong field; i.e., when frequency anticrossing induced by an interaction between the matter’s transition and a vacuum EM field becomes comparable to, or larger than, the transition and photon frequencies, the so-called ultrastrong coupling (USC) regime arises Forn-Díaz et al. (2018); Kockum et al. (2018).

This article is concerned with electron systems that are ultrastrongly coupled with photons in a cavity (or bosons with a resonant frequency, in more general, including magnons and phonons in long-wavelength limit). There are a variety of theoretically predicted many-body cavity quantum electrodynamics (QED) effects, states, and phases in the USC regime, which also offer possibilities for constructing unique protocols for ultrafast gates and ultrasecure state preparation for quantum information processing Forn-Díaz et al. (2018); Kockum et al. (2018). One of the most exciting aspects of cavity QED physics in the USC regime is the fact that the “light field” that the matter strongly couples with is not an external laser field but the vacuum cavity field, whose amplitude’s expectation value is zero. This fact distinctly differentiates USC physics from ordinary nonlinear optical phenomena, which are induced by a strong external field and thus inevitably involve excited and/or nonequilibrium states of matter. On the contrary, USC phenomena exclusively reflect the properties of the ground state of the matter–vacuum-field hybridized system in equilibrium. This new ground state, represented by a matter–vacuum-field entangled wavefunction Ciuti et al. (2005); Ashhab and Nori (2010); Felicetti et al. (2015), has characteristics that neither the original matter ground state (which can be metallic or insulating) nor the usual vacuum (which consists of fluctuating EM fields) possesses.

Specifically, in this article, we describe two systems shown in Fig. 1: (a) an ultrahigh-mobility two-dimensional electron gas (2DEG) in GaAs quantum wells placed inside a high-quality-factor($Q$) THz cavity Zhang et al. (2016); Li et al. (2018a) and (b) an ensemble of $\mathrm{Er}^{3+}$ spins interacting with $\mathrm{Fe}^{3+}$ magnons in a $\mathrm{ErFeO_{3}}$ crystal Li et al. (2018b). In both systems, we have observed USC — between the 2DEG and cavity photons in the former, and between electron spins and magnons in the latter. In order to describe the novel strong-field phenomena we observed in these systems, in Section II below, we present various types of Hamiltonians with different degrees of approximations appropriate for describing light–matter interactions in different coupling strength regimes. In the subsequent sections, Sections III and IV, we present more detailed models and results on the two systems. Finally, we summarize our findings and conclusions in Section V.

II Light–Matter Interaction Hamiltonians

In the presence of nonperturbative, or ultrastrong, coupling between light and matter, in general, we can no longer use the rotating-wave approximation (RWA) for describing light-matter interaction phenomena. Further, we can no longer neglect the ponderomotive energy $U_{\text{p}}$ Bucksbaum (1990); Chin et al. (2000); Srivastava et al. (2004); Wegener (2005). Even when the light is quite weak or even in the vacuum level, these kinds of approximations also fail when materials ultrastrongly interact with the EM fields Ciuti et al. (2005); Forn-Díaz et al. (2018); Kockum et al. (2018), i.e., when the interaction strength $g$ (which is equal to one half of the vacuum Rabi splitting) becomes comparable to, or larger than, the photon frequency $\omega_{a}$ and the transition frequency $\omega_{b}$ of the material.

In the case of $g\ll\omega_{a},\omega_{b}$, a system consisting of a photon mode with a resonance frequency $\omega_{a}$ and an ensemble of $N$ two-level atoms with a transition frequency $\omega_{b}$ can be described by the following Hamiltonian:

[TABLE]

In this paper, we call this the Tavis–Cummings Hamiltonian, following Ref. 6. Here, $\hat{a}$ ($\hat{a}^{\dagger}$) is the annihilation (creation) operator of a photon. $\hat{\sigma}_{i,z}=|\uparrow_{i}\rangle\langle\uparrow_{i}|-|\downarrow_{i}\rangle\langle\downarrow_{i}|$ is one of the Pauli operators representing the state of the $i$-th two-level atom, and $\hat{\sigma}_{i}=|\downarrow_{i}\rangle\langle\uparrow_{i}|$ is its lowering operator from the excited state $|\uparrow_{i}\rangle$ to the ground state $|\downarrow_{i}\rangle$. The last two terms in the Hamiltonian, $\hat{\sigma}^{\dagger}\hat{a}$ and $\hat{a}^{\dagger}\hat{\sigma}_{i}$, represent, respectively, the excitation of the atom through absorption of a photon and the relaxation of the atom through emission of a photon. $g/\sqrt{N}$ is the interaction strength per atom, while $g$ is that for the ensemble (collective excitation) of the atoms. In other words, the interaction strength $g$ scales as $\sqrt{N}$ with the increase in the number $N$ of atoms (precisely speaking, $g$ is proportional to the density of atoms beyond the long-wavelength approximation).

As discussed by R. H. Dicke in 1954 Dicke (1954), the above Hamiltonian is usually derived from the minimal-coupling Hamiltonian (in the Coulomb gauge) Cohen-Tannoudji et al. (1989):

[TABLE]

Here, $\hat{\bm{E}}_{\perp}(\bm{r})=-\hat{\bm{\varPi}}(\bm{r})/\varepsilon_{0}$ and $\hat{\bm{B}}(\bm{r})=\bm{\nabla}\times\hat{\bm{A}}(\bm{r})$ are the transverse electric field and magnetic flux density, respectively, expressed by the vector potential $\hat{\bm{A}}(\bm{r})$ and its conjugate momentum $\hat{\bm{\varPi}}(\bm{r})$ satisfying $[\hat{\bm{A}}(\bm{r}),\hat{\bm{\varPi}}(\bm{r}^{\prime})]=\mathrm{i}\hbar\bm{\delta}_{\perp}(\bm{r}-\bm{r}^{\prime})$, where $\bm{\delta}_{\perp}(\bm{r})$ is the transverse dyadic delta function Cohen-Tannoudji et al. (1989). $\hat{\bm{r}}_{i}$ and $\hat{\bm{p}}_{i}$ are the position and momentum of the $i$-th electron with charge $-e$ and mass $m$, satisfying $[\hat{\bm{r}}_{i},\hat{\bm{p}}_{j}]=\mathrm{i}\hbar\delta_{i,j}\bm{1}$. The last term $V(\{\hat{\bm{r}}\})$ represents the electrons’ one-body potential and many-body interactions in general. Here, we assume that, in this one-body potential, the oscillator strength is concentrated on a transition between the atomic ground state and an excited state with energy difference $\hbar\omega_{b}$. We also assume that the atoms are well isolated from each other and $V(\{\hat{\bm{r}}\})$ does not give many-body interactions. Further, we assume that the photons are confined in a cavity (embedding the atoms) with a resonance frequency $\omega_{a}$, and we focus only on one of the transverse directions of the polarization. Under these assumptions, in the long-wavelength approximation, $\hat{\mathcal{H}}_{\text{min}}$ can be rewritten as Dicke (1954); Rzażewski et al. (1975); Bamba and Ogawa (2014)

[TABLE]

In this paper, we call this the Dicke Hamiltonian. The last term is called the $A^{2}$, diamagnetic, or quadratic term. Its coefficient is usually derived as $D=g^{2}/\omega_{b}$ under the above assumptions. Here, $g^{2}$ is proportional to the oscillator strength or $d_{12}{}^{2}$ of the two-level atom, the atom density, and the vacuum fluctuation $\langle\hat{A}\hat{A}\rangle$ of the vector potential (in contrast to the ponderomotive energy Bucksbaum (1990); Chin et al. (2000); Srivastava et al. (2004); Wegener (2005), which is proportional to the intensity of light).

In the limit of $g\ll\omega_{a},\omega_{b}$, we can use the RWA (eliminating the counter-rotating terms $\hat{a}^{\dagger}\hat{\sigma}^{\dagger}_{i}$ and $\hat{a}\hat{\sigma}_{i}$) and neglect the last term in Eq. (3). Therefore, in that limit, $\hat{\mathcal{H}}_{\text{Dicke}}\approx\hat{\mathcal{H}}_{\text{Tavis--Cummings}}$. However, in the USC regime, where $g\gtrsim\omega_{a},\omega_{b}$, these approximations cannot be justified. Due to the presence of the counter-rotating terms, the expectation value $\langle\hat{a}^{\dagger}\hat{a}\rangle$ of the photon number becomes nonzero even in the ground state of the total system, and these photons are called virtual photons Ciuti et al. (2005). If the coefficient of the $A^{2}$ term satisfies $D<g^{2}/\omega_{b}$, we get a thermal second-order phase transition called the superradiant phase transition (SRPT) Hepp and Lieb (1973); Wang and Hioe (1973); Nataf and Ciuti (2010), where the expectation value $\langle\hat{a}\rangle$ of the EM field spontaneously becomes nonzero (not temporally oscillating) below a critical temperature in thermal equilibrium. In particular, when $D=0$, the SRPT occurs if the interaction strength satisfies $g^{2}>\omega_{a}\omega_{b}/4$ (i.e., the USC regime). However, considering the sum rule of the oscillator strengths, we usually get $D\geq g^{2}/\omega_{b}$, which forbids the SRPT Rzażewski et al. (1975); FN1 . A more general no-go theorem of the SRPT was also proposed in the minimal-coupling Hamiltonian, Eq. (2) Bialynicki-Birula and Rzażewski (1979); Gawȩdzki and Rza¸źewski (1981). Counter examples against that no-go theorem are still being discussed even in the current literature Grießer et al. (2016); Mazza and Georges (2019). However, the thermal SRPT has not yet been realized experimentally, while its nonequilibrium analogue has been demonstrated in cold atoms Baumann et al. (2010), and a thermal-equilibrium analogue has been theoretically proposed in an artificial system composed of a superconducting circuit Bamba et al. (2016).

In most experiments Forn-Díaz et al. (2018), the interaction strength $g$ and the coefficient $D$ of the $A^{2}$ term can be estimated through linear optical responses (absorption, transmission, or reflection spectra). In such cases, we can bosonize the ensemble of the atoms (through the lowest-level Holstein–Primakoff transformation), i.e., $\sum_{i=1}^{N}(\hat{\sigma}_{i,z}/2)\to\hat{b}^{\dagger}\hat{b}+\text{const.}$ and $\sum_{i=1}^{N}\hat{\sigma}_{i}/\sqrt{N}\to\hat{b}$, where $\hat{b}$ ($\hat{b}^{\dagger}$) is the annihilation (creation) operator of the collective excitation of the atoms. Then, the Dicke Hamiltonian is approximated as

[TABLE]

In this paper, we call this the Hopfield Hamiltonian. This approximation is justified in the weak excitation limit. The virtual photons appear also in the ground state of this Hamiltonian due to the counter-rotating terms. The SRPT cannot be described by the Hopfield Hamiltonian. Instead, we find an instability of the normal ground state (showing $\langle\hat{a}\rangle=\langle\hat{\sigma}_{i}\rangle=0$) under the same condition as that for the SRPT ($g^{2}>\omega_{a}\omega_{b}/4$ for $D=0$) FN2 .

For experimentally discussing the virtual photons and the existence of the SRPT, it is necessary to quantitatively evaluate the contributions of the counter-rotating terms and the $A^{2}$ term. In one of our experiments Li et al. (2018a), these contributions were evaluated separately by using a high-$Q$ THz cavity integrated with an ultrahigh-mobility 2DEG. As shown in Fig. 1(a), the THz cavity had a pair of silicon Bragg reflectors, each of which consisted of alternating layers of silicon and vacuum, and a defect thick silicon layer in the middle. The 2DEG, made of modulation-doped GaAs quantum wells grown by molecular-beam epitaxy, was attached to one of the surfaces of the defect layer; its position also overlapped with the position where the electric field reaches maximum at resonance. Under an external static magnetic flux density $B_{\text{stat}}$, the cyclotron motion of the 2DEG is excited resonantly by a THz wave with corotating circular polarization, which showed an anticrossing with the cyclotron resonance (CR) as we will see in Fig. 2. On the other hand, we found a frequency shift of the cavity mode with a circular polarization counter-rotating against the cyclotron motion (see the shaded area in Fig. 2). This shift is interpreted as the vacuum counterpart of the Bloch–Siegert (BS) shift Bloch and Siegert (1940), which usually appears due to an electron’s counter-rotating coupling with a strong light field. This vacuum BS shift was observed even when the average photon number was close to zero owing to the ultrastrong coupling, high-$Q$ cavity, and high-mobility 2DEG in our system. In order to describe the behavior of the resonance frequencies depending on the circular polarization, we need to consider explicitly the polarization degree of freedom of photons coupled with the CR, where the time reversal symmetry is broken under the static magnetic flux density $B_{\text{stat}}$. In contrast to the Hopfield Hamiltonian, Eq. (4), consisting of a single photon mode, the Hamiltonian in our system is expressed as

[TABLE]

Here, $\hat{a}_{+}$ ($\hat{a}_{-}$) is the annihilation operator of a photon with a circular polarization corotating (counter-rotating) against the cyclotron motion of the 2DEG in $B_{\text{stat}}$. $\hat{b}$ is the annihilation operator of the collective excitation of the 2DEG between the Landau levels with a transition frequency $\omega_{b}=e|B_{\text{stat}}|/m$. $\hat{a}_{+}$ interacts with $\hat{b}$ in the corotating manner as $\mathrm{i}\hbar g(\hat{b}^{\dagger}\hat{a}_{+}-\hat{a}^{\dagger}_{+}\hat{b})$, while $\hat{a}_{-}$ interacts with $\hat{b}$ in the counter-rotating manner as $\mathrm{i}\hbar g(\hat{b}^{\dagger}\hat{a}^{\dagger}_{-}-\hat{a}_{-}\hat{b})$. Because of this difference, we could clearly and quantitatively differentiate between the contributions of the counter-rotating terms (vacuum BS shift) and the $A^{2}$ terms, which will be discussed in Sec. III in addition to the derivation of $\hat{\mathcal{H}}_{\text{Landau}}$.

While the presence of a thermal SRPT in the minimal-coupling Hamiltonian (derived from the Maxwell equations and Newton’s equation of charged particles feeling the Lorentz force) in Eq. (2) is still under debate Bialynicki-Birula and Rzażewski (1979); Gawȩdzki and Rza¸źewski (1981); Grießer et al. (2016); Mazza and Georges (2019), there remains a possibility for realizing a thermal SRPT in systems where spins, instead of charges, interact with the EM fields, as pointed out by Knight et al. in 1978 Knight et al. (1978). A hint of such a possibility lies in certain magnetic phase transitions, which are caused by interactions between two species of spins within the same material. In particular, we focus on magnetic materials that can be modeled by the Dicke Hamiltonian under the replacement of photons with magnons in one of the spin species. If the magnetic phase transitions in such systems can be interpreted as magnon analogues of the SRPT, they will provide clues toward realizing the original photon SRPT. Recently, we have experimentally obtained evidence for ultrastrong and cooperative interactions between an ensemble of $\mathrm{Er}^{3+}$ ions and a magnon mode of the $\mathrm{Fe}^{3+}$ lattice in $\mathrm{ErFeO_{3}}$ crystals. As shown in Fig. 1(b), ordered $\mathrm{Fe}^{3+}$ spins form a magnon wave, while the electron paramagnetic resonance (EPR) of the $\mathrm{Er}^{3+}$ ions is interpreted as an atomic transition interacting with the magnon wave. Eventually, we can reduce a spin model of $\mathrm{ErFeO_{3}}$ into the Dicke Hamiltonian (with $D=0$) as

[TABLE]

Here, $\hat{a}$ is the annihilation operator of a magnon in the quasiferromagnetic (qFM) mode in the $\mathrm{Fe}^{3+}$ lattice. $\hat{\sigma}_{i}$ is the lowering operator between the lowest two $\mathrm{Er}^{3+}$ levels in a static magnetic field. As we will see in Fig. 4, we experimentally found an anticrossing between the $\mathrm{Er}^{3+}$ EPR and the qFM magnon mode, and the vacuum Rabi splitting ($2g$) reached a considerable fraction of the EPR and magnon frequencies. Further, by replacing $\mathrm{Er}^{3+}$ with $\mathrm{Y}^{3+}$ through chemical doping, we also found that the vacuum Rabi splitting $2g$ was proportional to the square root of the $\mathrm{Er}^{3+}$ density. These are discussed in Sec. IV in addition to the derivation of the Dicke Hamiltonian in Eq. (6) from the spin model of $\mathrm{ErFeO_{3}}$.

III Cyclotron resonance interacting with THz photons

As shown in Fig. 1(a), we consider a planar cavity embedding a 2DEG parallel to the cavity mirrors. Generalizing the minimal-coupling Hamiltonian in Eq. (2), we explicitly consider the cavity mirror structure described by a position-dependent relative dielectric constant $\varepsilon_{\text{cav}}(z)$, where $z$ is the direction perpendicular to the mirrors and the 2DEG. For simplicity, we focus only on the EM wave without an in-plane wavenumber. The total Hamiltonian is expressed as

[TABLE]

Following the quantization of the EM wave in an inhomogeneous dielectric medium in Ref. 32, we describe the Hamiltonian of the EM wave (excluding the external static field) as

[TABLE]

The conjugate momentum of the vector potential $\hat{A}_{\xi}(z)$ corresponds to the electric displacement field in the cavity medium (but without the 2DEG) as $\hat{\varPi}_{\xi}(z)=-\varepsilon_{0}\varepsilon_{\text{cav}}(z)\hat{E}_{\xi}(z)$. These operators satisfy $[\hat{A}_{\xi}(z),\hat{\varPi}_{\xi^{\prime}}(z^{\prime})]=\mathrm{i}\hbar\delta_{\xi,\xi^{\prime}}\delta(z-z^{\prime})/S$, where $S$ is the area of the $x$–$y$ plane. On the other hand, the kinetic energy of the 2DEG is expressed as

[TABLE]

Here, $\hat{\bm{\pi}}_{i}\equiv\hat{\bm{p}}_{i}+e\bm{A}_{\text{stat}}$ and $\bm{A}_{\text{stat}}$ gives the external static magnetic flux density $\bm{B}_{\text{stat}}=\bm{\nabla}\times\bm{A}_{\text{stat}}$. $z_{\text{2DEG}}$ is the position of the 2DEG inside the cavity. $N=n_{\text{2DEG}}S$ is the number of electrons for surface density $n_{\text{2DEG}}$ of the 2DEG. We do not consider electron-electron Coulomb interactions, since they do not affect the linear optical response of the 2DEG due to Kohn’s theorem Kohn (1961).

Let us first rewrite the EM Hamiltonian $\hat{\mathcal{H}}_{\text{EM}}$ in terms of the annihilation and creation operators of the photons. From the Heisenberg or Hamilton’s equations (giving the Maxwell equations) derived from $\hat{\mathcal{H}}_{\text{EM}}$, we get a wave equation for the vector potential $\hat{A}_{\xi}(z)$ as

[TABLE]

Following the quantization procedure in Ref. 32, by replacing $\hat{A}_{\xi}(z,t)$ with $f_{k}(z,t)/\sqrt{\varepsilon_{\text{cav}}(z)}$ and by performing the temporal Fourier transformation, we rewrite the wave equation as

[TABLE]

Solving this equation, we can determine the eigenfrequencies $\{\omega_{k}\}$ and eigenfunctions $\{f_{k}(z)\}$ of the EM wave in the cavity. $k$ is the mode index. The wavefunctions are normalized as $\int\mathrm{d}z\ f_{k}(z)f_{k^{\prime}}(z)=\delta_{k,k^{\prime}}$. They also satisfy the completeness as $\sum_{k}f_{k}(z)f_{k}(z^{\prime})=\delta(z-z^{\prime})$. From the complete set of these eigenmodes, the operators of the vector potential and its conjugate momentum are described as

[TABLE]

Here, $\hat{a}_{k,\xi}$ is the annihilation operator of a photon in the $k$-th mode with a linear polarization in the $\xi=x,y$ direction, satisfying $[\hat{a}_{k,\xi},\hat{a}^{\dagger}_{k^{\prime},\xi^{\prime}}]=\delta_{k,k^{\prime}}\delta_{\xi,\xi^{\prime}}$. The Hamiltonian of the EM wave is then rewritten as

[TABLE]

In general, there are not only the cavity modes (spatially localized modes) but also continuous modes (transmission modes).

On the other hand, the 2DEG shows the CR with a frequency $\omega_{\text{cyc}}=e|B_{\text{stat}}|/m$. Introducing the lowering operator $\hat{c}_{i}\equiv(\hat{\pi}_{i,y}+\mathrm{i}\hat{\pi}_{i,x})/\sqrt{2m\hbar\omega_{\text{cyc}}}$ between the Landau levels satisfying $[\hat{c}_{i},\hat{c}^{\dagger}_{j}]=\delta_{i,j}$ Elvang and Polchinski (2003), we rewrite the kinetic energy of the 2DEG as

[TABLE]

The last term is the $A^{2}$ term. The second term contains the lowering and raising processes involved with the non-Hermitian vector potential

[TABLE]

where the annihilation operator of the $\pm$ circularly polarized photon in the $k$-th mode is defined as

[TABLE]

Introducing the bosonic operator of a collective excitation of the 2DEG coherently interacting with the EM wave

[TABLE]

we can rewrite the interaction term as

[TABLE]

where the interaction strength for mode $k$ is expressed as

[TABLE]

On the other hand, the $A^{2}$ term is rewritten as

[TABLE]

Then, when we focus only on the CR collective excitation described by $\hat{b}$ and one cavity mode with resonance frequency $\omega_{a}=\omega_{k}$, the total Hamiltonian is finally expressed as Eq. (5).

Figures 2(a) and (b) show numerically calculated THz transmission spectra for different values of the static magnetic flux density $B_{\text{stat}}$. The interaction strengths are assumed to be (a) $g/2\pi=37.5\;\mathrm{GHz}$ and (b) $g/2\pi=75.0\;\mathrm{GHz}$ (see the details in the caption). The dashed red and black lines show the resonance frequencies of the cavity mode and CR, respectively, without considering their interaction. Figure 2(c) shows experimental data, where we estimated $g/2\pi=150.1\;\mathrm{GHz}$. We used a circularly polarized THz wave as a probe beam, and it was sufficiently weak so that any optical nonlinear effects did not appear. At $B_{\text{stat}}>0$, the CR and the THz wave are corotating, and an anticrossing (vacuum Rabi splitting) appears. In contrast, at $B_{\text{stat}}<0$, the CR is counter-rotating against the THz wave, and a frequency shift is obtained, instead of the anticrossing, as highlighted with the gray shaded areas.

These anticrossing and frequency shift can be reproduced by the Hamiltonian $\hat{\mathcal{H}}_{\text{Landau}}$ in Eq. (5). In its third term, the CR interacts with the $+$ and $-$ circularly polarized cavity modes in the corotating manner $\mathrm{i}\hbar g(\hat{b}^{\dagger}\hat{a}_{+}-\hat{a}^{\dagger}_{+}\hat{b})$ and the counter-rotating manner $\mathrm{i}\hbar g(\hat{b}^{\dagger}\hat{a}^{\dagger}_{-}-\hat{a}_{-}\hat{b})$, respectively. This aspect is clearly reflected in the eigenfrequencies $\omega_{\pm}$ of the coupled modes (polariton modes) with the $\pm$ circular polarization determined by

[TABLE]

Here, $\tilde{g}\equiv g\sqrt{\omega_{a}/\omega_{\text{cyc}}}$. These equations are obtained from the Fourier transform of the Heisenberg equations derived from $\hat{\mathcal{H}}_{\text{Laudau}}$. The second term in Eq. (22b) corresponds to the optical susceptibility. Since it is proportional to $(\omega_{+}-\omega_{\text{cyc}})^{-1}$, the $+$ circularly polarized photon resonantly interacts with the CR, and the anticrossing is obtained as seen in Fig. 2 for $B_{\text{stat}}>0$. On the other hand, the $-$ circularly polarized mode does not show such an anticrossing since the optical susceptibility is proportional to $(\omega_{-}+\omega_{\text{cyc}})^{-1}$. Since this frequency dependence purely comes from the counter-rotating coupling between CR and $-$ circularly polarized photons, the frequency shift observed for $B_{\text{stat}}<0$ in Fig. 2 (gray shaded areas) can be interpreted as a BS shift. In contrast to the standard BS shift proportional to the intensity of the EM wave, the frequency shift in our system depends on the vacuum fluctuation $\langle\hat{A}_{\xi}(z_{\text{2DEG}})\hat{A}_{\xi}(z_{\text{2DEG}})\rangle=\hbar f(z_{\text{2DEG}})^{2}/[2\varepsilon_{0}\varepsilon_{\text{cav}}(z_{\text{2DEG}})\omega_{a}S]=\hbar\tilde{g}^{2}m/(Ne^{2}\omega_{a})$ of the EM field, and it is interpreted as the vacuum counterpart of the BS shift.

In this way, the contributions of the corotating and counter-rotating couplings can be distinguished through the circular polarization in the CR system. Such a distinguishability cannot be obtained for the Lorentz-type susceptibility

[TABLE]

which is derived from the Hopfield Hamiltonian, Eq. (4), for $D=g^{2}/\omega_{b}$ and $\tilde{g}=g\sqrt{\omega_{a}/\omega_{b}}$. In this case, the corotating contribution $(\omega_{b}-\omega)^{-1}$ and the counter-rotating one $(\omega_{b}+\omega)^{-1}$ cannot be clearly separated in the linear optical spectra.

On the other hand, the second term in Eq. (22a) comes from the $A^{2}$ term. In the case of $B_{\text{stat}}=0$, we get $\omega_{\text{cyc}}=g=0$ while keeping $\tilde{g}\neq 0$. Then, the last term in Eq. (22a) disappears, and we get $\omega_{\pm}(B_{\text{stat}}=0)=\sqrt{\omega_{a}{}^{2}+2\tilde{g}^{2}}$, whose frequency shift $\Delta\omega\equiv\omega_{\pm}(B_{\text{stat}}=0)-\omega_{a}$ from the bare cavity frequency $\omega_{a}$ is purely the contribution of the $A^{2}$ term. In the limit of $B_{\text{stat}}\to-\infty$, this frequency shift is canceled by the vacuum BS shift as seen in Fig. 2. In this way, the contribution of the $A^{2}$ term is also clearly distinguished from those of the corotating and counter-rotating couplings in the CR system.

The requirement for observing the vacuum Rabi splitting is that $2g$ should be larger than the line broadening. Owing to this relatively easy requirement, the ultrastrong vacuum Rabi splitting has been observed in a variety of systems Forn-Díaz et al. (2018); Kockum et al. (2018). In contrast, the maximum vacuum BS shift is $\Delta\omega=\sqrt{\omega_{a}{}^{2}+2\tilde{g}^{2}}-\omega_{a}$, and it is approximately expressed as $\Delta\omega\approx\tilde{g}^{2}/\omega_{a}$ in the case of $\tilde{g}\ll\omega_{a}$. Therefore, in order to observe the vacuum BS shift, in addition to the ultrastrong $g$ in the CR system and the circularly polarized THz wave, the high-$Q$ THz cavity and the high-mobility 2EDG are essential for realizing small enough line broadening as shown in Fig. 2(c). These conditions were all satisfied in our experiments Li et al. (2018a).

IV Cooperative interaction in magnetic materials

In this section, we will derive the Dicke Hamiltonian, Eq. (6), from the following spin model for $\mathrm{ErFeO_{3}}$:

[TABLE]

Following the discussion by G. F. Herrmann in 1963 Herrmann (1963), we express the Hamiltonian of the $\mathrm{Fe}^{3+}$ spins by a two-sublattice model as

[TABLE]

Here, $\hat{\bm{S}}_{i}^{A/B}$ is the operator of the $i$-th $\mathrm{Fe}^{3+}$ spin vector with magnitude $S=5/2$ in its $A/B$ sublattice. $\sum_{\mathrm{n.n.}}$ means a summation over all the nearest neighbor couplings. $N$ is the number of $\mathrm{Fe}^{3+}$ spins in each sublattice, i.e., there are in total $2N$ spins representing the $\mathrm{Fe}^{3+}$ subsystem. $J_{\text{Fe}}$ and $D_{\text{Fe}}$ are the isotropic and antisymmetric exchange constants, respectively. $A_{x}$, $A_{z}$, and $A_{xz}$ are the magnetic anisotropy energies. On the other hand, we assume that the Hamiltonian of the $\mathrm{Er}^{3+}$ spins is simply expressed under an external static magnetic flux density $\bm{B}_{\text{stat}}$ as

[TABLE]

Here, $\hat{\bm{\mu}}_{i}$ is the operator of the $i$-th $\mathrm{Er}^{3+}$ magnetic moment vector expressed by the Pauli operators $\{\hat{\sigma}_{i,\xi}\}$ as $\hat{\bm{\mu}}_{i}=-\mu_{\text{B}}(\mathfrak{g}_{x}\hat{\sigma}_{i,x},\mathfrak{g}_{y}\hat{\sigma}_{i,y},\mathfrak{g}_{z}\hat{\sigma}_{i,z})^{\mathrm{t}}$, where $\mathfrak{g}_{\xi}$ is the anisotropic $g$-factor. In the following, the $\mathrm{Er}^{3+}$ spin vector is expressed as $\hat{\bm{R}}_{i}=(\hat{\sigma}_{i,x},\hat{\sigma}_{i,y},\hat{\sigma}_{i,z})^{\text{t}}$. We do not consider the $\mathrm{Er}^{3+}$–$\mathrm{Er}^{3+}$ exchange interactions. Finally, the $\mathrm{Fe}^{3+}$–$\mathrm{Er}^{3+}$ exchange interactions are expressed as

[TABLE]

Here, $J_{A/B}$ and $\bm{D}_{A/B}$ are the symmetric and antisymmetric exchange constants, respectively. In this way, we assume that the $\mathrm{Fe}^{3+}$ spins in the $i$-th unit cell of $\mathrm{ErFeO_{3}}$ are represented by $\hat{\bm{S}}^{A/B}_{i}$ (two-sublattice model), each of which is in fact a sum of two real $\mathrm{Fe}^{3+}$ spins in the four-sublattice model. On the other hand, in the same manner as Ref. 12, the four $\mathrm{Er}^{3+}$ spins in the $i$-th unit cell are represented simply by one Pauli vector $\hat{\bm{R}}_{i}$.

In the ground state of the $\mathrm{Fe}^{3+}$ subsystem, their spins are ordered antiferromagnetically along the $z$ axis as depicted in Fig. 3, but slightly canted toward the $x$ direction by angle $\beta_{0}=(-1/2)\arctan[(A_{xz}+D_{\text{Fe}})/(J_{\text{Fe}}-A_{x}+A_{z})]$ Herrmann (1963). The equilibrium spins are expressed as $\bar{\bm{S}}^{A}_{i}=S(\sin\beta_{0},0,-\cos\beta_{0})^{\mathrm{t}}$ and $\bar{\bm{S}}^{B}_{i}=S(\sin\beta_{0},0,\cos\beta_{0})^{\mathrm{t}}$. We define the modulations $\{\delta\hat{S}_{\ell,T},\delta\hat{S}_{\ell,Y}\}$ from them as depicted in Fig. 3, and these operators satisfy $[\delta\hat{S}_{\ell,T},\delta\hat{S}_{\ell^{\prime},Y}]=\mathrm{i}\delta_{\ell,\ell^{\prime}}$. The spin modulations are then expressed as

[TABLE]

Extending Herrmann’s calculation Herrmann (1963) into a magnon model propagating in the $z$ direction (under averaging in the $x$–$y$ plane) Tsang et al. (1978), we can derive the equations of motion for these modulations as

[TABLE]

Here, $\gamma=\mathfrak{g}\mu_{\text{B}}/\hbar$ is the gyromagnetic ratio for the free electron $g$-factor $\mathfrak{g}$ and the Bohr magneton $\mu_{\text{B}}$. The coefficients are expressed as Herrmann (1963)

[TABLE]

where $z=6$ is the number of neighboring $\mathrm{Fe}^{3+}$ sites for each $\mathrm{Fe}^{3+}$ spin. Then, the Hamiltonian of the $\mathrm{Fe}^{3+}$ spins is approximated (bosonized) as

[TABLE]

Here, $N_{z}$ and $2N_{z}$ are the number of unit cells and of $\mathrm{Fe}^{3+}$ spins, respectively, in the $z$ direction. In terms of the annihilation operator $\hat{a}_{K}$ of a magnon with a dimensionless wavenumber $K$, satisfying $[\hat{a}_{K},\hat{a}^{\dagger}_{K^{\prime}}]=\delta_{K,K^{\prime}}$, the modulation operators are expressed as

[TABLE]

The Hamiltonian in Eq. (31) is rewritten as

[TABLE]

where the eigenfrequency is expressed as

[TABLE]

Here, $K=0$ and $K=\pi$ correspond to the qFM and quasi-antiferromagnetic (qAFM) modes Herrmann (1963), respectively.

Dividing the $\mathrm{Fe}^{3+}$ spin operators into the equilibrium value $\bar{\bm{S}}^{A/B}$ and the modulation operators $\{\delta\hat{\bm{S}}^{A/B}_{i}\}$, the rest of the Hamiltonian is rewritten as

[TABLE]

where

[TABLE]

We also averaged the $\mathrm{Er}^{3+}$ spins in the $x$–$y$ plane. In the experiment of Ref. 12, the main contribution to the EPR of $\mathrm{Er}^{3+}$ ions was the the external static magnetic field along the $c$ ($z$) axis, and the exchange contribution from the equilibrium $\mathrm{Fe}^{3+}$ spins $\bar{\bm{S}}^{A/B}$ was a minor one. Then, the EPR Hamiltonian is approximated as

[TABLE]

where $\omega_{\text{EPR}}=2g_{z}\mu_{\text{B}}B_{\text{stat}}/\hbar$. A more detailed analysis was performed in the Supplementary Materials of Ref. 12. On the other hand, the EPR–magnon interaction is dominated by the symmetric exchange interactions, since $J_{A/B}\gg|\bm{D}_{A/B}|$ is usually satisfied. Due to a symmetry analysis, we can find $J_{A}=J_{B}=J$ Li et al. (2018b). Then, when we assume the long-wavelength limit and focus only on the qFM magnon mode with frequency $\omega_{\text{FM}}=\omega_{K=0}$, which shows the anticrossing with EPR as seen in Fig. 4(a), the interaction Hamiltonian is approximated as

[TABLE]

Since the last term gives a nonlinear optical response for the approximated EPR Hamiltonian in Eq. (38), the interaction Hamiltonian is approximated in the linear responses as

[TABLE]

In experiments, we can dilute the density of $\mathrm{Er}^{3+}$ ions on demand by substitutional doping with nonmagnetic $\mathrm{Y^{3+}}$ ions. Considering the dilution of $\mathrm{Er}^{3+}$ spins by factor $x$ due to the replacement with $\mathrm{Y}^{3+}$ ions and also the thermal distribution of electron spins in $\mathrm{Er}^{3+}$ ions, the net fraction $\eta_{\text{spin}}\equiv x\tanh(-\hbar\omega_{\text{EPR}}/k_{\text{B}}T)$ of the EPR-contributing $\mathrm{Er}^{3+}$ density was introduced in Eq. (40) Li et al. (2018b). The interaction strength can be expressed as

[TABLE]

In this way, the interaction strength $g$ is proportional to the square root of $\eta_{\text{spin}}$. From Eqs. (33), (38), and (40), we get the Dicke Hamiltonian, i.e., Eq. (6). Note that no $A^{2}$ term appears in this Dicke Hamiltonian (even without the renormalization mentioned in FN1 , since the magnon–EPR interaction is derived from the $\mathrm{Er}^{3+}$–$\mathrm{Fe}^{3+}$ exchange interactions in the spin model of Eq. (24), not from the kinetic energy in the minimal-coupling Hamiltonian.

Figure 4(a) shows experimentally observed absorption peak frequencies of $\mathrm{ErFeO_{3}}$ at 10 K as a function of static magnetic flux density $B_{\text{stat}}$ along the $c$ ($z$) axis. The dashed red and black lines show uncoupled EPR and magnon frequencies, respectively. The interaction strength $g$ can be directly read from the vacuum Rabi splitting $2g$ in the graph, and we found $g/\omega_{\text{FM}}=14.7\;\%$, reaching the ultrastrong coupling regime. Moreover, by measuring multiple samples of $\mathrm{Er}_{x}\mathrm{Y}_{1-x}\mathrm{FeO_{3}}$ at multiple temperatures, we found that $g$ showed a proportionality with $(\eta_{\text{spin}}\omega_{\text{FM}})^{1/2}$ as shown in Fig. 4(b) and as theoretically derived in Eq. (41). This proportionality is experimental evidence of the Dicke cooperativity between the ensemble of $\mathrm{Er}^{3+}$ ions and the magnon mode in the $\mathrm{Fe}^{3+}$ lattice. While the origin of their interaction is short-range (nearest-neighboring) exchange interactions as modeled in Eq. (27), the EPR of the $\mathrm{Er}^{3+}$ ensemble cooperatively interacts with the $\mathrm{Fe}^{3+}$ magnon mode, which propagates through the $\mathrm{Fe}^{3+}$–$\mathrm{Fe}^{3+}$ exchange interactions.

While $\omega_{\text{FM}}$ were slightly modified depending on the sample and temperature, we assume $J_{\text{Fe}}=4.96\;\mathrm{meV}$ as reported in Ref. 37. From the proportionality $2g/(2\pi)\approx(\eta_{\text{spin}}\omega_{\text{FM}}/2\pi)\times 0.238\;\mathrm{THz}^{1/2}$ obtained in Fig. 4(b), we can estimate the $\mathrm{Fe}^{3+}$–$\mathrm{Er}^{3+}$ symmetric exchange interaction strength as $J=2.95\;\mathrm{meV}$.

V Summary

The counter-rotating light–matter coupling and the $A^{2}$ (quadratic) term are essential for exploring the virtual photons and the SRPT in the USC regime. In the CR system, apart from the contribution of the corotating coupling (vacuum Rabi splitting), those of the counter-rotating coupling and $A^{2}$ term can be clearly distinguished, respectively, as the vacuum BS shift and the cavity frequency shift at zero static magnetic field. In a bulk $\mathrm{ErFeO_{3}}$ crystal, the coupled system of the $\mathrm{Er}^{3+}$ spins and $\mathrm{Fe}^{3+}$ magnons can be described by the Dicke Hamiltonian. If we find magnon analogues of the SRPT in such magnetic materials, it would give us a hint of realizing the original photon SRPT. The possibility of SRPT still remains in the light–spin coupling systems experimentally, while the thermal SRPT has been discussed mostly in ideal theoretical models. The systems without the time reversal symmetry and a variety of magnetic materials would advance the exploration of virtual photons and the SRPT in the USC regime.

ACKNOWLEDGMENTS

This research was supported by the National Science Foundation (Cooperative Agreement DMR-1720595), the U.S. Army Research Office (grant W911NF-17-1-0259), the JST PRESTO program (grant JPMJPR1767), KAKENHI (grant 26287087), and the ImPACT Program of the Council for Science, Technology and Innovation (Cabinet Office, government of Japan). We thank Katsumasa Yoshioka for his help with all measurements done with circularly polarized THz radiation. Saeed Fallahi, Geoff Gardner, and Michael Manfra provided the 2DEG sample for the Landau polariton work in Section III. Ning Yuan, Maolin Xiang, Kai Xu, and Shixun Cao provided the $\mathrm{Er}_{x}\mathrm{Y}_{1-x}\mathrm{FeO}_{3}$ crystals for the work described in Section IV.