Ultrastrong magnon–magnon coupling dominated by antiresonant interactions

Exotic quantum vacuum phenomena are predicted in cavity quantum electrodynamics systems with ultrastrong light-matter interactions. Their ground states are predicted to be vacuum squeezed states with suppressed quantum fluctuations owing to antiresonant terms in the Hamiltonian. However, such predictions have not been realized because antiresonant interactions are typically negligible compared to resonant interactions in light-matter systems. Here we report an unusual, ultrastrongly coupled matter-matter system of magnons that is analytically described by a unique Hamiltonian in which the relative importance of resonant and antiresonant interactions can be easily tuned and the latter can be made vastly dominant. We found a regime where vacuum Bloch-Siegert shifts, the hallmark of antiresonant interactions, greatly exceed analogous frequency shifts from resonant interactions. Further, we theoretically explored the system’s ground state and calculated up to 5.9 dB of quantum fluctuation suppression. These observations demonstrate that magnonic systems provide an ideal platform for exploring exotic quantum vacuum phenomena predicted in ultrastrongly coupled light-matter systems.


Supplementary Note 1: Rice Advanced Magnet with Broadband Optics (RAMBO)
RAMBO is schematically illustrated in Supplementary Fig. 1. Magnetic fields are generated using a 30 T pulsed mini-coil with a 12 mm inner diameter. The magnet is housed in a custom-built liquid nitrogen cryostat that serves to prevent Joule heating from changing the resistance of the coil between magnetic field pulses and to lower the resistance of the coil, thereby increasing the peak magnetic field strength. Coaxial electrodes connect the mini-coil to a 5.6 mF capacitor bank capable of being charged up to over 1, 800 V, as is required for 30 T magnetic field pulses. Single crystals of YFeO 3 are held in the center of the mini-coil by being mounted on the narrow end of a tapered sapphire pipe (minimum inner diameter 6 mm, minimum outer diameter 8 mm). The sapphire pipe is attached to the copper coldfinger of a liquid helium cryostat by an indium ring.
The coldfinger is covered by a radiation shield to prevent heating from the cryostat walls. With these measures in place, the lowest sample temperature achievable is roughly 12 K, although the measurements presented in this work were done without liquid helium and were therefore done at room temperature. The liquid helium cryostat and liquid nitrogen cryostat are combined by a bellows joint, enabling fine adjustments of the sample's position in the bore of the magnet.
Optical access to the sample is provided through the tapered sapphire pipe and windows on the liquid helium and liquid nitrogen cryostats. The tapering of the sapphire pipe increases the optical access to the sample. The magnetic field profile is determined by measuring the current delivered to the mini-coil using a Rogowski Coil (Powertek, CWT 60 Mini Rogowski Coil). The relationship between the current delivered to the mini-coil and the magnetic field profile was calibrated using a pickup coil placed at center of the bore of the magnet. To determine the temporal relationship between laser pulses from our amplified Ti:Sapphire and the magnetic field profile, a photodiode is placed adjacent to the THz generation crystal to measure scattered light. By using single-shot THz detection, we are able to measure transmitted THz electric fields at three different magnetic field strengths for a single magnetic field pulse, as illustrated in Fig. 1c of the main text.

Supplementary Note 2: Single-shot Detection
The transmitted THz waveform is obtained from the image of the reflective echelon following our procedure in Reference 1 . For each column of pixels imaged on the CMOS camera, the total intensity of one polarization component is calculated by vertically integrating over one half of the column, and the total intensity of the other polarization component is calculated by vertically integrating over the other half of the column. The difference of these two intensities is labeled ∆I. Although the optical probe pulse is magnified before the echelon to ensure uniform intensity over the image, we normalize both polarization components' intensities by the total intensity of the column (I o ) to normalize for nonuniform probe intensity over the echelon. The ratio ∆I/I o , is proportional to the THz electric field strength. This process is repeated for all columns in order to find the THz electric field strength at each horizontal pixel location on the camera, as indicated in We start from the treatment given by Herrmann 2 by considering interactions between the two Fe 3+ spin sublattices. The free energy of the system, normalized by the magnitude M 0 of the two spin sublattices S 1 and S 2 , was first derived by Herrmann for the Γ 4 phase with an applied magnetic field along the c-axis. The normalized free energy accounting for a tilted magnetic field in the b − c plane is given by: where axis. Note that we assume the anisotropy constant A xz to be zero for simplicity, as in References 3-5 .
In equilibrium, in the Γ 4 phase, R 1 and R 2 are ordered antiferromagnetically along the a-axis, with a slight canting towards the c-axis. In terms of the ferromagnetic vector F = R 1 + R 2 and the antiferromagnetic vector G = R 1 − R 2 , we have that F points along the c-axis, and G points along the a-axis.
The equations of motion for R i are given by: where γ is the gyromagetic ratio and ∇ i V is the gradient of V with respect to (X i , Y i , Z i ). The equations of motion for R i take a simpler form when transformed from Cartesian coordinates (1, 0, 0) in equilibrium. The transformation, as illustrated in Supplementary Fig. 4, is given by: The equilibrium positions are found by solving the equations of motion forṘ i = 0.
where ... represents an evaluation of the derivative at the equilibrium position. For notational simplicity, we rewrite Supplementary Eq. 5 as: which defines the constants a, b, c, d, e, f , g, and h. These are given by: a = E(cos 2 β z cos 2β y − sin 2 β z ) + D sin 2β z cos β y + 2A xx cos 2 β z cos 2 β y Given that δT i and δY i are canonical variables, the equations of motion given in Supplementary Eq. 5 are recovered from Hamilton's equations: for the following Hamiltonian: The equations of motion and their physical interpretation are more illuminating when rewrit-ten in terms of F and G. The transformation is given by: and the transformed Hamiltonian is given by: where the coefficients are functions of the spin sublattice equilibrium positions and derivatives of the free energy, and are given by: Since we get from Supplementary Eq. 10, Supplementary Eq. 11, [δT i , δY j ] = iδ i,j , and [δT i , δT j ] = [δY i , δY j ] = 0, we consider δF x and δF y /2 sin β z (δG x and δG y /2 sin β z ) to be canonically conjugate variables.
Hamilton's equations for δF x and δF y /2sinβ z are given by: and similarly for δG x and δG y /2sinβ z : which can be written more compactly as: These equations of motion yield two magnonic eigenfrequencies: as discussed in Methods.
To fully understand magnonic interactions, we proceed to quantize our Hamiltonian. Here, we introduce magnonic creation and annihilation operators for each magnon mode:â,â † for the generalized qFM mode andb,b † for the generalized qAFM mode. We may rewrite fluctuations in F in terms of these new operators: Similarly, we may rewrite fluctuations in G in terms ofb andb † : We can rewrite our Hamiltonian, Supplementary Eq. 12, in terms of the magnonic creation and annihilation operators: where ω 0a and ω 0b are generalized qFM and qAFM magnon frequencies given by: and g 1 and g 2 are the co-rotating and counter-rotating coupling strengths, respectively, given by: We confirm our quantized Hamiltonian by calculating its predicted magnon frequencies and comparing it to our previously derived magnon frequencies. The magnon frequencies are obtained from the equations of motion for the creation and annihilation operators: which yields the following equations of motion: Expressions for the two positive magnonic eigenfrequencies are given by:

Supplementary Note 4: Hopfield-Bogoliubov Transformation
Here, we diagonalize our Hamiltonian using a Hopfield-Bogoliubov transformation. We introduce coupled magnon annihilation operatorsB L (B U ) describing the LM (UM), which are expressed in terms of the generalized qFM (qAFM) operatorsâ (b) by: for j = L, U . The coefficients are constrained by the relation: and are determined by the following eigenvalue problem:  where Ω + (Ω − ) corresponds to the UM (LM) eigenfrequency and is provided in Supplementary Eq. 28.

Supplementary Note 5: Evaluating Squeezing
The variances ofXĉ ,φ andXd ,φ , defined in Methods, can be easily evaluated in the coupled magnon ground state |0 by inverting the Hopfield-Bogoliubov transformation given in Supplementary Eq. 29 and rewriting the quadratures in terms ofB j . The fluctuation inXĉ ,φ in the coupled magnon ground state is given by: and a similar expression can be derived for 0|(∆Xd ,φ ) 2 |0 .
Supplementary Note 6: Origin of g 1,2 and relation to magnetic parameters Here, we discuss the origin of the magnon-magnon coupling, g 1,2 , and their relation to magnetic parameters, E, D, A xx , A zz .
First, we clarify that the origin of g 1,2 is the tilted magnetic field. As discussed in the main text, when H DC is zero or applied along the c-axis, the two spin sublattices maintain π rotational symmetry about the c-axis. Given fixed and opposite parities of the qFM and qAFM modes under this symmetry, their hybridization is prevented in this geometry. When β y is nonzero, their hybridization is no longer forbidden, where β y defines how far the spins are pulled from the a − c plane ( Supplementary Fig. 4). This is also clear from the analytical definitions of the coupling strengths ( Supplementary Eqs. 24, 25), where it follows from Supplementary Eq. 7 and Supplementary Eq. 13 that g 1,2 exactly vanishes when β y is zero. Therefore, nonzero β y is the origin of our magnon-magnon coupling, and we achieve this condition by applying a strong magnetic field with a nonzero component along the b-axis.
Second, we discuss the relation between magnetic parameters and the anisotropy of g 1 and g 2 . Large magnitudes of g 2 are crucial for the exotic quantum vacuum phenomena predicted in the ultrastrong coupling regime, and thus there is great motivation for understanding how the strength of these counter-rotating interactions depend on magnetic parameters. Supplementary Fig. 12 shows numerical calculations for normalized coupling strengths |g 1,2 |/ω 0 vs. θ, where ω 0 is the frequency at which the generalized qFM and qAFM modes cross (also defined in Fig. 3c of the main text). In each plot, we calculate |g 1,2 |/ω 0 vs. θ as we tune one magnetic parameter. Specifically, in each plot the dotted lines correspond to scaling the chosen magnetic parameter by 2, the dashed lines correspond to scaling the chosen magnetic parameter by 0.5, and the solid lines correspond to the true value of the chosen magnetic parameter (i.e. scaled by 1). In Supplementary   Fig. 12a, 12b, 12c, and 12d, we tune E, A xx , D, and A zz , respectively.
From Supplementary Fig. 12, we see that g 2 /g 1 becomes larger as E, A xx increase and as D, A zz decrease i.e. the coupling becomes more anisotropic. Similarly, we see that g 2 /g 1 becomes smaller as E, A xx decrease and as D, A zz increase i.e. the coupling becomes more isotropic. The relationship between the anisotropicity and the magnetic parameters becomes more illuminating when we recall that increasing E, A xx and decreasing D, A zz reduces the spin canting angle (β z in Supplementary Fig. 4), whereas decreasing E, A xx and increasing D, A zz increases the spin canting angle. The relationship between the spin canting angle and magnetic parameters can be analytically derived by solving for the steady states of Supplementary Eq. 2. When H DC is applied along the b-axis, one can show that the expression is given by: where H y is the b-component of H DC . Therefore, the coupling becomes more isotropic as the spins become more canted, and the coupling becomes more anisotropic as the spins become more antiparallel. Finally, we highlight that the origin of our ground state squeezing is the counterrotating terms. Therefore, as the coupling becomes more anisotropic due to stronger counterrotating interactions the degree of squeezing will also be amplified.   Supplementary Fig. 13. Coupling strength variation with θ. a, Plots of |g 1 |/2π vs. θ for different fixed magnetic fields, and similarly for |g 2 |/2π vs. θ in b. Importantly, we observe that |g 2 |/2π is a monotonically increasing function of θ.