Quantum Simulation of an Extended Dicke Model with a Magnetic Solid

The Dicke model describes the cooperative interaction of an ensemble of two-level atoms with a single-mode photonic field and exhibits a quantum phase transition as a function of light--matter coupling strength. Extending this model by incorporating short-range atom--atom interactions makes the problem intractable but is expected to produce new phases. Here, we simulate such an extended Dicke model using a crystal of ErFeO$_3$, where the role of atoms (photons) is played by Er$^{3+}$ spins (Fe$^{3+}$ magnons). Through magnetocaloric effect and terahertz magnetospectroscopy measurements, we demonstrated the existence of a novel atomically ordered phase in addition to the superradiant and normal phases that are expected from the standard Dicke model. Further, we elucidated the nature of the phase boundaries in the temperature--magnetic-field phase diagram, identifying both first-order and second-order phase transitions. These results lay the foundation for studying multiatomic quantum optics models using well-characterized many-body condensed matter systems.

The Dicke model describes the cooperative interaction of an ensemble of two-level atoms with a single-mode photonic field and exhibits a quantum phase transition as a function of lightmatter coupling strength.Extending this model by incorporating short-range atom-atom interactions makes the problem intractable but is expected to produce new physical phenomena and phases.Here, we simulate such an extended Dicke model using a crystal of ErFeO 3 , where the role of atoms (photons) is played by Er 3+ spins (Fe 3+ magnons).Through terahertz spectroscopy and magnetocaloric effect measurements as a function of temperature and magnetic field, we demonstrated the existence of a novel atomically ordered phase in addition to the superradiant and normal phases that are expected from the standard Dicke model.Further, we elucidated the nature of the phase boundaries in the temperature-magnetic-field phase diagram, identifying both first-order and second-order phase transitions.These results lay the foundation for studying multiatomic quantum optics models using well-characterized many-body solid-state systems.

Main text
The Dicke model in quantum optics describes the cooperative, coherent coupling of an ensemble of two-level atoms with a single-mode light field 1 .Despite its simplicity, the model hosts a rich variety of phenomena that are significant in diverse contexts, such as cavity quantum electrodynamics 2 , condensed matter physics 3 , and quantum information science 4,5 .A prominent feature of the model is a second-order quantum phase transition (QPT), known as the superradiant phase transition (SRPT), which occurs when the light-matter coupling strength, g, exceeds a threshold 6,7 .
When the system enters the superradiant phase, atomic and photonic polarizations spontaneously emerge, producing a unique many-body ground state that enables studies of unusual light-matter entanglement 8 , two-mode squeezed states [9][10][11] , and quantum chaos 12 .
Although the atomic ensemble in the original Dicke model was assumed to be noninteracting, it has been known from the early days that atom-atom interactions are important for explaining, for example, the dephasing and intensity correlation functions of fluorescent spectra 13,14 .Hence, there has long been interest in extending the Dicke model to include an atom-atom interaction (represented by strength J); see Fig. 1.Such an extended Dicke model, or the g-J model, should display an interplay of two types of interatomic interactions -i.e., the photonic-field-mediated long-range interaction, and the direct short-range interaction.Intuitively, one can expect the ground state of the system to crucially depend on the ratio g/J, with a superradiant phase (an atomically ordered phase) favored for large (small) g/J.However, no analytical solutions can be obtained for the g-J model, motivating one to simulate it using a well-characterized many-body quantum system.
Computational studies of the g-J model under various approximations have revealed an array of new phenomena, such as a first-order QPT [15][16][17][18][19] , a shift of the SRPT boundary 20,21 , amplification of the integrablity-to-chaos transition 22 , modifications of matter-matter entanglement 20,23 , and alteration of the nature of an excited-state QPT 18,24 .To examine these phenomena, several experimental platforms, including atomic Bose-Einstein condensates 25,26 , superconducting qubits 27,28 , and quantum dots 15 , have been proposed as quantum simulators, but successful simulations have not been achieved.
Here, we present a novel protocol of using a crystal of erbium orthoferrite (ErFeO 3 ), an antiferromagnetic (AFM) insulator, as a solid-state quantum simulator of the g-J model.The magnetic properties of ErFeO 3 are governed by the moments carried by the Er 3+ and Fe 3+ spin subsystems and their interplay 29 .A previous study has revealed Dicke cooperativity in the Er 3+ -Fe 3+ interaction 30 , demonstrating the resemblance of the magnetic Hamiltonian of ErFeO 3 to the Dicke Hamiltonian.Namely, the paramagnetic Er 3+ ions (the magnons of ordered Fe 3+ spins) play the role of the atomic ensemble (light field), and the spin-magnon interaction is formally similar to the g-term in the Dicke model.What further strengthens this analogy is a magnetic phase transition of the crystal that exhibits all traits that would be expected for a Dicke SRPT.When the temperature (T ) becomes lower than 4 K, the Er 3+ lattice develops C-type AFM order 31 (with the ferromagnetic chains along z), and a zone-boundary Fe 3+ magnon mode condenses, displacing the staggered moments away from the x-z plane 32,33 ; this corresponds to the emergence of atomic and photon polarizations in the standard SRPT.In Bertaut's notation, the magnetic transition is of the Γ 2 → Γ 12 type (Fig. 2a).Mean-field calculations using a realistic spin model captures the simultaneous order parameter (OP) onset of both the Er 3+ and Fe 3+ spin components, ⟨Σ − z ⟩ and ⟨S y ⟩, respectively (Fig. 2b), indicating that the Γ 2 → Γ 12 transition is a magnonic SRPT 34 , with the Γ 2 and Γ 12 phases corresponding to the normal (N) and superradiant (S) phases, respectively.
One way to observe the OP onset is to monitor the quasi-antiferromagnetic (qAFM) magnon mode of Fe 3+ spins through terahertz (THz) time-domain spectroscopy, which has been utilized to reveal the configuration of Fe 3+ ions in rare-earth orthoferrites 35 .By performing THz transmission measurements on a z-cut ErFeO 3 crystal in the Faraday geometry, we obtained absorption coefficient (α) spectra, derived from the imaginary part of the refractive index (Supplementary Section 2), as a function of T , as shown in Fig. 2c.The observed bright absorption line is the qAFM mode, which has been thoroughly studied in previous studies 35 .It is the evolution of this mode in distinct phases of the g-J model that is of interest throughout this study.A continuous OP-like onset, or a kink, is observed at the N → S transition boundary (< 4 K, blue dashed line).The frequency shift of the qAFM magnon mode in the S phase from that in the N phase is thus a sensitive reporter of the qAFM magnon condensate density, namely, the Fe 3+ OP of the S phase.
The J-term is inherently built into the magnetic Hamiltonian of ErFeO 3 since the Er 3+ -Er 3+ exchange interaction, albeit being weak, is known to be present 36 .Spectroscopic measurements have also revealed a fine frequency splitting within the Er 3+ electron paramagnetic resonance lines 30 , which is attributable to the Er 3+ -Er 3+ exchange interaction.The presence of both the gand J-terms sets the stage for ErFeO 3 to simulate the g-J model.Nonetheless, although the g-term-driven S phase can find correspondence to the Γ 12 phase in ErFeO 3 , the g/J ratio set for the crystal stipulates that a pure atomic (A) phase, which is driven exclusively by the J-term, would not appear in equilibrium.For ErFeO 3 , the A phase would be an Er 3+ ordered phase without involving any OP onset in the Fe 3+ subsystem.Therefore, to achieve quantum simulation of the g-J model, we must search for a way to invoke an explicit A phase through an S → A transition.
Our theoretical consideration suggests that subjecting ErFeO 3 to a static magnetic field (H) along the z axis can potentially induce an S → A transition.This can be understood by writing the simplified magnetic Hamiltonian (Supplementary Section 1) in the second-quantized form as where a two-sublattice approximation is adopted for both Er 3+ and Fe 3+ for a total of N 0 unit cells.Here, ω π , â † π , and âπ are the energy, creation and annihilation operators for the Fe 3+ qAFM magnon mode, respectively; ω Er is the frequency of Er 3+ spins as two-level systems at H = 0; , where g z is the Landè g factor, µ B is the Bohr magneton, and µ 0 is the vacuum permeability, is the H-induced Zeeman frequency of Er 3+ ; and g and J are the Er 3+ -magnon and Er 3+ -Er 3+ coupling strengths, leading to the gand J-terms of the g-J Hamiltonian, respectively.Σp = 2N 0 i=1 σi,p /2, where σp are Pauli matrices and p ∈ {x, y, z}, is the collective Er 3+ spin operator, with its superscript "+" ("−") denoting the sum (difference) of the two sublattices.
The way these operators appear in Eq. ( 1) is crucial for interpreting the ground-state energetics.
Specifically, the g-term features a product of the Fe 3+ magnon field operator i(â † π − âπ ) and the Σ− z component of Er 3+ spins, thereby favoring antiparallel alignment of Er 3+ sublattices and Fe 3+ magnon condensation in the S phase (the onsets of ⟨ Σ− z ⟩ and ⟨S y ⟩ in Fig. 2b), whereas the J-term This would tip the balance between the g-term and the J-term, since Σ− x appears only in the J-term but not in the g-term.
As shown in Fig. 3a, an S → A transition is indeed recovered in the calculated mean-field phase diagram of the spin Hamiltonian (Supplementary Section 1) within the T -H parameter space, for T < 2.8 K, with a critical field ranging from 0.35 T to 0.5 T, depending on T .Increasing the field to above 1 T and elevating T to above 4 K would both push the system across the thermodynamic phase boundary into the N phase.A triple point (at 2.8 K and 0.5 T, decorated by a yellow star) marks the location where the S, A, and N phases converge.Figure 3b shows the calculated normalized spin components as the OPs of the magnetic phases, for a line cut along the H axis at T = 0 K, traversing sequentially the S → A and the A → N boundaries.We identify that the Fe 3+ OP, represented by ⟨S y ⟩, is finite in the S phase but near-zero in the A phase.The Er 3+ OP, on the other hand, is finite in both the S and A phases, but undergoes a switch from the Further, the OP evolution indicates that the S → A boundary is an abrupt-type, first-order phase transition, while the A → N boundary is a continuous-type, second-order phase transition.
Summarizing the mean-field calculation results, Fig. 3c pictorially shows the predicted Fe 3+ and Er 3+ spin order in each phase.Starting from the N phase, the two sublattices of Fe 3+ are antiparallel along z with zero y-component, while Er 3+ spins remain paramagnetic (no order).The A phase is characterized by Fe 3+ order that is identical to that of the N phase, but the Er 3+ subsystem develops canted AFM order where the sublattice moments are antiparallel along x (⟨ Σ− x ⟩ ̸ = 0), with canting along z (⟨ Σ+ z ⟩ ̸ = 0).In the S phase, the Er 3+ order takes the ⟨ Σ+ configuration, and the staggered moment of the Fe 3+ sublattices undergoes a rotation about the x axis, bringing its y-component to nonzero.
The S → A transition can be considered as a spin-flop transition in terms of Er 3+ ions.One conventional way to characterize the transition is to monitor the magnetic susceptibility through which the existence of the AFM ordering of Er 3+ spins in the A phase has been previously observed 37 , although the configuration of Fe 3+ spins was left ambiguous.Our magnetization measurements showed clear S → A and S → N phase boundaries (squares in Fig. 4a and Supplementary Section 2).However, a strong and non-uniform demagnetizing effect that broadens the phase boundary 38 likely prevented us from clearly identifying the A → N phase boundary.This is because the shape of our sample for magnetization measurements was a thin irregularly shaped disk cut from the sample used for THz measurements, rather than a sphere, which would have produced a uniform demagnetizing field inside the sample 37 .Nonetheless, a disk-shaped sample with a large lateral size was necessary for performing THz transmission measurements.
To demonstrate the A → N phase boundary, i.e., the breakdown of the AFM order of Er 3+ spins, we performed magnetocaloric effect (MCE) experiments which are sensitive to the magnetic entropy landscape of a material.Namely, the Grüneisen ratio 39 measures the slope of isentropes in the T − H plane 40 .Since the heat capacity C H is always a positive quantity, the sign of ( ∂T ∂H ) S is always opposite to ( ∂S ∂H ) T .Furthermore, sharp changes in entropy S due to phase transitions will appear as step functions in ( ∂T ∂H ) S 40 , or peaks if is measured in a quasi-adiabatic environment 41 .Thus, by measuring the differential change in sample temperature with respect to the magnetic field, ( ∂T ∂H ) S , the T − H magnetic phase diagram can be measured.We note, that the demagnetization factor can have a small effect on the MCE measurements 42 , namely that temperature and field shifts can occur, but the qualitative features should be present.
The T -H phase diagram of ErFeO 3 , and the obtained results are summarized in Fig. 4a.We configured a MCE measurement in a Physical Property Measurement System in the quasi-adiabatic condition 43 , and took raw data traces of sample temperature variation versus magnetic field at a ramping rate of 5 × 10 −3 T/s with dH > 0 (Supplementary Section 2); the sensitivity of temperature variation of our instrument reached 5 × 10 −4 K. To identify H-induced phase transitions, the first-order derivative ( ∂T ∂H ) S was approximated as dT /d(µ 0 H) (Supplementary Section 2), whose local extremes correspond to the transition boundaries 41 .The traces clearly exhibit two maxima for T < 2.8 K, corresponding to the S→A (red dashed line in Fig. 4a) and A → N (blue dashed line) boundaries, and one maximum for 2.8 K < T < 4 K, corresponding to the S → N (blue dashed line) boundary.These results are qualitatively consistent with the T − H phase diagram reported previously 37 , where quantitative shifts likely come from demagnetization effects.
Once we experimentally investigated the evolution of the atomic ensemble (or the Er T -dependence of α spectra at select µ 0 H values, respectively.We found that the bright absorption lines can be assigned to either the quasi-ferromagnetic (qFM) mode or the qAFM mode 35 , the latter of which can be an OP for the Fe 3+ spins.
In the H-dependent color map at 1.4 K (Fig. 4b), three lines are observed.The lowest frequency line, which does not pick up intensity until 0.8 T, is the qFM mode, while the other two are both qAFM magnons, albeit belonging to distinct phases.The middle (upper) line, which is located at 0.8 THz at 0 T (1 THz at 0.5 T), is the qAFM mode of the S (A & N) phase.The S → A transition can be identified to occur at 0.5 T (red dashed line), where the upper line emerges.The qAFM magnons belonging to the S and A phases coexist within 0.5 T < µ 0 H < 1 T, consistent with the prediction that the S → A transition is of first order and is thus inhomogeneous, until the middle line vanishes at >1 T (blue dashed line) owing to entrance into the N phase.The 3.2 K map (Fig. 4c) shows a different behavior; the qAFM magnon (0.88 THz at 0 T) of the S phase continuously shifts to connect with that of the N phase in frequency, forming an OP-like onset for µ 0 H < 0.7 T (blue dashed line), signaling a second-order N → S transition boundary.Such a frequency shift is absent in the 4.4 K map (Fig. 4d) since the N phase persists throughout the whole T -dependent color maps at constant H further corroborate our assignments of the phase transitions.Again from the 0 T map (Fig. 4e), a continuous OP-like onset of the qAFM mode shift is observed (< 4 K, blue dashed line).This echoes Fig. 4c in showing the continuous nature of the N → S transition, and establishes that the frequency shift of the qAFM magnon in the S phase from that in the N phase can be the Fe 3+ OP of the S phase.Intriguingly, this OP is demonstrated to be zero in the A phase.We read this fact from the 0.75 T map (Fig. 4f), for which an N → A transition is expected upon lowering T .Although a residual mode pertaining to the S phase exists (as mentioned earlier when discussing Fig. 4b), the qAFM mode (unlabeled line) frequency does not undergo any noticeable OP-like anomaly across the N → A transition; it is as featureless as the qAFM mode within the 1.25 T map (Fig. 4g), for which the N phase persists throughout the whole T range.This unambiguously demonstrates that the spin order in the A phase only involves Er 3+ ordering but not any Fe 3+ OP, consistent with our expectation depicted in Fig. 3. Finally, phase boundaries determined by the THz experiments are overlaid (as solid circles) on top of the MCE phase diagram in Fig. 4a, showing overall agreement.
A potential impact of this analogy is the possibility of being applied to other members of the rare-earth orthoferrite family or orthochromite compounds.For example, spin-reorientation phase transitions 44,45 (Γ 4 → Γ 2 ) in RFeO 3 (R=Yb, Er, and Tb) would mimic the SRPT.In YbFeO 3 , where the Yb 3+ -Yb 3+ interaction (J) is negligible, it would be a potential playground for studying the standard Dicke model (g model).At the boundaries of the phase transition of YbFeO 3 , the qFM mode of Fe 3+ shows a kink, and a transition inside the ground doublet of Yb 3+ ions shows a softening 46 .This simultaneous kink and softening is one of the hallmarks of a magnonic SRPT 47 .
It was also suggested that TbFeO 3 can be regarded as the magnetic phase transition of the Jahn-Teller type 48,49 that would resemble a magnonic SRPT.In ErFeO The advantages of using the low-temperature phase transition, as opposed to the ∼80 K spinreorientation phase transition, of ErFeO 3 in simulating the extended Dicke model can be summa-rized as follows.First, the low-temperature phase transition allows us to simulate the first-order phase transition into the A phase, which is the main point of this work and does not exist in the spinreorientation transition at 87 K. Second, since we deal with the lowest two energy levels (Kramers doublet) of Er 3+ ions in the low-temperature phase transition, theoretical analysis is directly relevant to the Dicke model, compared to the multiple crystal-field energy levels involved in the 87 K phase transition.Third, and most importantly, at high temperatures, thermally populated magnons are not negligible.Such thermal magnons will prevent studies of the vacuum magnons responsible for the Dicke superradiant phase transition, which occurs in thermal equilibrium without any external driving.For example, one consequence of the superradiant phase transition induced by vacuum bosonic fields is a two-mode perfect squeezed vacuum at the critical point 34 .A finite number of thermally excited magnons will mask such interesting quantum phenomena.
In summary, through THz magnetospectroscopy and magnetocaloric effect experiments, we studied a crystal of ErFeO 3 to simulate the g-J model, which is an extended Dicke model that includes not only the bosonic-field-mediated long-range interatomic interactions but also direct short-range interactomic interactions.In addition to the superradiant and normal phases expected from the standard Dicke model, we identified a new phase, an atomic phase, which is driven by the short-range J-term in the Hamiltonian.Further, we elucidated the nature of the various phase boundaries, distinguishing between first-order and second-order transitions.These results demonstrated the potential of ErFeO 3 as a simulator of quantum optics Hamiltonians.More specifically, in the context of Dicke physics, this condensed matter platform may lead to the possibilities of assisting quantum chaos 22 and modifying matter-matter entanglement 20,23 with tunability given through an external magnetic field.Bridging the gap between quantum optics and many-body correlated physics, our results will find broader application in the design of hybrid quantum systems with superior controllability, such as the Dicke-Ising machine 27,50 and the Dicke-Lipkin-Meshkov-Glick model 18,23 .Furthermore, the ability to transition between the superradiant and atomic phases via a nonthermal knob provides opportunities to study unconventional quantum criticality 51 and chaos-assisted thermalization 52 .

Methods
Methods, including statements of data availability and any associated accession codes and references, are available in Supplementary Information.
3+ spins) in the extended Dicke Hamiltonian, we turned to elucidate the photonic counterpart (or the Fe 3+ spins).The ambiguity of the configuration of Fe 3+ spins and the nature of the transition boundaries require us to monitor the qAFM magnon mode of Fe 3+ spins in THz magnetospectroscopy experiments.Unlike the static measurements, responses from different domains in the A phase can be distinguished in the frequency domain, illuminating the nature of the phase transition.The measurements were performed within the same T -H parameter space as that of the MCE experiments.Figures 4b-d and Figs.4e-g show the H-dependence of α spectra at select T values and the

3 ,
where a crystal field transition (∼1.5 THz) is responsible for the spin-reorientation transition (T = 87 K), the crystal field levels would play the role of an ensemble of two-level atoms in the Dicke model.To prove that Dicke physics is at work, however, one must show Dicke cooperativity, i.e., the coupling strength g must exhibit cooperative enhancement g ∝ √ N , where N is the number of two-level atoms.In addition, mapping their spin Hamiltonians into the Dicke models is required to establish this analogy.No attempts have been made to develop an analogy between the spin-reorientation transition and the Dicke superradiant phase transition.

Figure 1 :Figure 2 :
Figure 1: The extended Dicke model, or the g-J model, where an ensemble of interacting two-level

Figure 4 :
Figure 4: Mapping out the T -H phase diagram of ErFeO 3 in H ∥ z. a, Phase boundaries determined by