Controlling the magnetic state of the proximate quantum spin liquid α-RuCl₃ with an optical cavity

Abstract

Harnessing the enhanced light-matter coupling and quantum vacuum fluctuations resulting from mode volume compression in optical cavities is a promising route towards functionalizing quantum materials and realizing exotic states of matter. Here, we extend cavity quantum electrodynamical materials engineering to correlated magnetic systems, by demonstrating that a Fabry-Pérot cavity can be used to control the magnetic state of the proximate quantum spin liquid α-RuCl₃. Depending on specific cavity properties such as the mode frequency, photon occupation, and strength of the light-matter coupling, any of the magnetic phases supported by the extended Kitaev model can be stabilized. In particular, in the THz regime, we show that the cavity vacuum fluctuations alone are sufficient to bring α-RuCl₃ from a zigzag antiferromagnetic to a ferromagnetic state. By external pumping of the cavity in the few photon limit, it is further possible to push the system into the antiferromagnetic Kitaev quantum spin liquid state.

Introduction

The realization of magnetic van der Waals (vdW) materials with thicknesses down to the monolayer limit has sparked a new interest in fundamental aspects of two-dimensional magnetism^1,2,3. Due to a competition of strong anisotropy, fluctuations, and spin-orbit effects, vdW materials are prime candidates to host exotic phenomena such as topological phase transitions, magnetic skyrmions, and quantum spin liquids^4,5. In addition, the electronic, magnetic, and optical properties of these materials are sensitive to a wide range of material engineering techniques such as strain^6,7, nanostructuring⁸, electric fields^9,10, and moiré twisting^11,12, allowing their state to be tuned with high precision.

Recent progress has also established optical engineering techniques as a method to functionalize quantum materials and to reach exotic (out-of-equilibrium) topological phases^{13,14,15,16,17}. However, driving a system with lasers is associated with excessive heating when the frequency becomes multi-photon resonant with electronic transitions^18,19. A way to circumvent this problem is to embed the system in an optical cavity, where the effective light-matter coupling is enhanced via mode volume compression and the state of the material can be modified in an equilibrium setting^{20,21,22,23,24,25,26,27}. Due to the strong interaction between light and charged excitations, polaritonic control of material and chemical properties via the cavity vacuum has so far focused on electronic and phonon mediated phase transitions^28,29. Alternatively, to address magnetic systems, optical cavities have been combined with external drives^30,31, or been used to modify a system’s excited state properties^32,33. While currently efforts are made to extend the cavity framework to a broader class of materials, and to construct a unified first principles description of cavity quantum fluctuations and quantum matter^20,23,34, experiments demonstrating polaritonic control of materials are scarce^35,36. Therefore, to transform this promising approach into a powerful experimental tool, it is of key importance to identify candidate materials where cavity engineering techniques can be explored.

Here, we extend the concept of cavity quantum electrodynamics (c-QED) engineering into the magnetic regime and identify such a candidate system, by demonstrating how an optical cavity can be used to control the magnetic ground state of the proximate quantum spin liquid α-RuCl₃ via shaping the quantum fluctuations of the cavity^20,23,29,37. Depending on the cavity frequency, photon occupation, and the strength of the effective light-matter coupling, we find that it is possible to transform the equilibrium zigzag antiferromagnetic order into any of the magnetic phases supported by the extended Kitaev model (see Eq. (1) and Fig. 2). As a key result we find that for frequencies of a few THz and for moderate light-matter couplings, the interaction between the magnetic system and the vacuum fluctuations of the cavity is sufficient to transform α-RuCl₃ from a zigzag antiferromagnet to a ferromagnet. In contrast to the meta-stable states obtained by driving the system with classical light, the magnetic state resulting from the interaction with the quantum fluctuations of the cavity is a true equilibrium state denoted the photo ground state (PGS)²⁴. Pumping the cavity in the few photon regime, it is further possible to push the system into the Kitaev quantum spin liquid state and to retrieve the non-equilibrium phase diagram of the semi-classical limit³⁸. Our results pave the way for utilizing c-QED to induce and control long-lived exotic states in quantum materials.

Results

Low-energy model and equilibrium magnetic order

The vdW material α-RuCl₃ consists of layers of Ru atoms arranged in an hexagonal lattice and surrounded by edge-sharing octahedra of Cl ions (Fig. 1a). Due to the crystal field the Ru d-orbitals are split by an energy Δ_cf into a lower e_g and a higher t_2g manifold, and the strong spin-orbit coupling further splits the t_2g states by an energy Δ_soc into a lower j_eff = 3/2 quartet, fully occupied in the ground state, and a j_eff = 1/2 doublet with a single hole. The Cl ions are assumed to be completely filled in the ground state, with their p-orbitals separated from the Ru d-orbitals by a charge-transfer energy Δ_pd. The Ru energy level structure is schematically shown in Fig. 1c. The local interactions of the t_2g manifold are described by a Hubbard-Kanamori Hamiltonian H_U^38,39, which takes into account the intra-orbital Hubbard interaction U, the inter-orbital interaction ${U}^{{\prime} }$, the Hund’s coupling J_H and the spin-orbit coupling λ (see “Methods” for a discussion of the model). The Ru holes can either hop directly between Ru atoms or move via the Cl ligands, as described by the kinetic Hamiltonian H_t. In the strong coupling limit t/U ≪ 1 virtual hopping processes give rise to effective magnetic interactions via the superexchange mechanism. Due to the orbital alignment a Kitaev interaction arises from ligand-mediated hopping over 90^∘ bond angles⁴⁰, while sub-dominant exchange and anisotropy terms arise from direct Ru–Ru interactions.

Fig. 1: Crystal and magnetic structure of monolayer α-RuCl₃.

a Crystal structure of monolayer α-RuCl₃, with the magnetic Ru ions (orange) arranged in a hexagonal lattice and exhibiting a zigzag antiferromagnetic order. The surrounding octahedra of Cl ions give rise to a crystal field splitting of the Ru d-orbitals (see panel c) b The magnetic system interacts with a strength g with a cavity electric field of frequency Ω. c Energy level structure of the magnetic Ru ions, leading to an effective j_eff = 1/2 magnetic moment in the Ru t_2g manifold. d Magnetic interactions of Ru moments, colored according to the Kitaev bonds x, y, and z.

Experimentally α-RuCl₃ is found to exhibit a zigzag antiferromagnetic order below the Néel temperature T_N ≈ 7 K, as indicated in Fig. 1. At zero temperature, first principles calculations show that this zigzag state is approximately degenerate with a ferromagnetic state⁴¹, and might only be stabilized by spin quantum fluctuations⁴². In addition, signatures of a Kitaev quantum spin liquid (QSL) state have been found upon applying an external magnetic field along the out-of-plane direction^43,44. Together these results indicate that although α-RuCl₃ orders at low temperatures, its ground state is proximate to several competing magnetic orders and the magnetic phase diagram is determined by a delicate competition of different magnetic interactions. This makes α-RuCl₃ an excellent candidate material to explore the competition between cavity and spin quantum fluctuations.

In the following the material will be assumed to have C₃ symmetry, which is satisfied to a very good degree in α-RuCl₃. All parameters of the local and kinetic Hamiltonians H_U and H_t were calculated from first principles as discussed in the “Methods” section, and give a zigzag antiferromagnetic ground state in line with observations.

Light-matter coupling in a cavity

Within the low-energy description discussed above, the main effect of the cavity is to modify the hopping amplitudes of the Ru holes. Inside the cavity the total Hamiltonian is given by $H={H}_{U}+{H}_{t}+{\sum }_{\lambda }\hslash {\Omega }_{\lambda }{\hat{n}}_{\lambda }$, where ℏΩ_λ is the energy of photon mode λ and ${\hat{n}}_{\lambda }$ is the corresponding number operator. The photon field modifies the kinetic Hamiltonian H_t by introducing the replacements ${\hat{c}}_{i\alpha \sigma }^{{\dagger} }{\hat{c}}_{i\beta \sigma }\to {e}^{i{\phi }_{ij}}{\hat{c}}_{i\alpha \sigma }^{{\dagger} }{\hat{c}}_{i\beta \sigma }$, where the Peierls phases ϕ_ij are proportional to the quantum vector potential $\hat{{{{\bf{A}}}}}={\sum }_{\lambda }({g}_{\lambda }{{{{\bf{e}}}}}_{\lambda }{\hat{a}}_{\lambda }^{{\dagger} }+{g}_{\lambda }^{* }{{{{\bf{e}}}}}_{\lambda }^{* }{\hat{a}}_{\lambda })$. For a given mode λ the bare light-matter coupling is defined as ${g}_{\lambda }=ea/\sqrt{2\epsilon \hslash {\Omega }_{\lambda }V}$, where e is the elementary charge, ℏ the reduced Planck constant, ϵ the relative permittivity and V the cavity mode volume. For a two-dimensional Fabry-Pérot cavity the lowest energy photon mode has a frequency Ω₀ = πc/L_z, and assuming a cavity with area A measured in in units of the squared Ru–Ru distance a², the bare light-matter coupling of this mode is ${g}_{0}=e/\sqrt{2\pi \epsilon \hslash cA}\approx 0.12/\sqrt{A}$.

Similarly we define an effective light-matter coupling $g=e/\sqrt{2\pi \epsilon \hslash c{A}_{{{{\rm{eff}}}}}}$ of a single cavity mode by introducing an effective mode area A_eff, which accounts for dielectric properties of the environment^{23,24,45,46,47}. In addition, we consider an effective single-mode approximation obtained by integrating over photon modes with momenta in a disk of radius q_m around q = 0^48,49. The light-matter coupling of the effective mode is proportional to the photonic density of states at q = 0, which grows as $\rho \sim \sqrt{A}$ in the macroscopic limit and cancels the $1/\sqrt{A}$ scaling of the bare coupling. This single effective mode will therefore have a coupling ${g}_{{{{\rm{eff}}}}}=\rho g\approx 0.12/\sqrt{{A}_{c,{{{\rm{eff}}}}}}$, where A_c,eff is the effective unit cell mode area. Although the value of A_c,eff can in principle be fixed by comparisons to experiment, since it is hard to estimate we will in the following vary g_eff to address the effect of different light-matter coupling strengths.

Effective coupled spin-photon model

The total Hamiltonian is down-folded to the spin sector by eliminating H_t to fourth order in virtual ligand-mediated processes using quasi-degenerate perturbation theory. The structure of the perturbation expansion allows for the down-folding to be performed separately within each given photon sector, resulting in a coupled spin-photon Hamiltonian of the form

$${{{\mathcal{H}}}}=\mathop{\sum}\limits_{{{{\bf{nm}}}}}\left({{{{\mathcal{H}}}}}_{s,{{{\bf{nm}}}}}+{\delta }_{{{{\bf{nm}}}}}\mathop{\sum}\limits_{\lambda }\hslash {\Omega }_{\lambda }{n}_{\lambda }\right)\left\vert {{{\bf{n}}}}\right\rangle \left\langle {{{\bf{m}}}}\right\vert .$$

(1)

Here ${{{{\mathcal{H}}}}}_{s,{{{\bf{n}}}}m}$ is the spin Hamiltonian in the sector connecting cavity number states $\left\vert {{{\bf{n}}}}\right\rangle$ and $\left\vert {{{\bf{m}}}}\right\rangle$, where m = {m₁, m₂, …, m_L} for L modes. The Hamiltonian in each photon sector is

$$\begin{array}{ll}{{{{\mathcal{H}}}}}_{s,{{{\bf{nm}}}}}\,=\,\mathop{\sum}\limits_{\langle ij\rangle }\left(\begin{array}{ccc}{S}_{i}^{\alpha }&{S}_{i}^{\beta }&{S}_{i}^{\gamma }\end{array}\right){\left(\begin{array}{ccc}J&\Gamma &{\Gamma }^{{\prime} }\\ \Gamma &J&{\Gamma }^{{\prime} }\\ {\Gamma }^{{\prime} }&{\Gamma }^{{\prime} }&J+K\end{array}\right)}_{{{{\bf{nm}}}}}\left(\begin{array}{c}{S}_{j}^{\alpha }\\ {S}_{j}^{\beta }\\ {S}_{j}^{\gamma }\end{array}\right)\\ \qquad\qquad\quad+\,{B}_{{{{\bf{nm}}}}}\mathop{\sum }\limits_{i}{\hat{{{{\bf{e}}}}}}_{B}\cdot {{{{\bf{S}}}}}_{i}\end{array}$$

(2)

where the magnetic interactions J, K, Γ, and ${\Gamma }^{{\prime} }$ all depend on the light-matter coupling g = {g₁, g₂, …, g_L} and the photon numbers n and m. Each bond 〈ij〉 is labeled by the indexes αβ(γ) ∈ {xy(z), yz(x), zx(y)} and is denoted a γ-bond in accordance with Fig. 1d. The induced magnetic field B points along the [111] direction of the local spin axes, as defined by the unit vector ${\hat{{{{\bf{e}}}}}}_{B}$.

The form of the spin-photon Hamiltonian in Eq. (2) is valid as long as the C₃ symmetry is preserved. Since a single linearly polarized mode breaks rotational symmetry, such a mode induces additional terms in the spin Hamiltonian. The rotational symmetry breaking is also reflected in a strong dependence of the cavity-modified spin parameters on the polarization direction, as shown in Supplementary Figure 1. However, in a Fabry-Pérot cavity the lowest mode is doubly degenerate with two possible in-plane polarizations, and including both these modes in the description restores the rotational symmetry of ${{{\mathcal{H}}}}$. The C₃ symmetry can also be maintained by considering a single circularly polarized mode in a chiral cavity²³, as will be done in the following.

We note that in principle the Hamiltonian also contains couplings to collective spin excitations, that could become important in the macroscopic limit and in presence of long range order. However, since the photon momentum at optical energies is small (q ≈ 0), such a coupling would mainly involve the ferromagnetic order parameter. In most of the phase space this coupling is strongly suppressed, and we therefore neglect such terms in the following. To quantitatively address their effects would require the development of new computational techniques, e.g., based on tensor network or quantum Monte Carlo methods.

Cavity dissipation

Below we discuss how the magnetic phase diagram of α-RuCl₃ is modified both by a dark and driven cavity, by assuming a perfect cavity limit. In any realistic cavity, photon losses will limit the applicability of our results, since after a typical time-scale κ⁻¹ (where κ is the cavity decay rate) a typical photon will have been emitted. For the results obtained for the dark cavity dissipation imposes no additional restriction, since the dark cavity is expected to be insensitive to cavity losses. To validate this, we have checked that the magnetic transitions in this regime arise mainly from changes of the magnetic parameters within the $\left\vert 0\right\rangle \left\langle 0\right\vert$ photon sector. For the results of the driven cavity, the light-induced change in the magnetic parameters is found to be of the order 1–10 meV, and thus a κ ≈ 1 meV or smaller would be necessary to observe the cavity-induced magnetic phase transitions. This corresponds to a cavity quality factor of Q = Ω/κ ≈ 1000 or higher in the optical frequency regime, which is achievable with present technologies.

Magnetic phases of the cavity photo ground state

The photo ground state (PGS) and magnetic phase diagram of the coupled spin-photon system was obtained by exact diagonalization of Eq. (1) on the 24-site spin cluster shown in Fig. 2a, including a single effective cavity mode with circular polarization. This spin cluster is the minimal one known to respect all sub-lattice symmetries of the magnetic system⁵⁰. The calculation results in a true equilibrium state where the cavity is close to its vacuum state and the magnetic system is in the state favored by the photo-induced interactions. To perform the calculations we assume a perfect cavity, while a more quantitative description should account for dissipation. However, as discussed above, the results below are expected to be robust against such effects.

Fig. 2: Magnetic phases of the photo ground state.

a Spin cluster employed for the exact diagonalization studies. b Paths through the magnetic phase diagram traced out by the system as the light-matter coupling g_eff increases from g_eff,i = 0 to g_eff,f = 0.5. The paths are colored according to the cavity photon energy ℏΩ, and the small arrows show the direction of the paths. The polar coordinate denotes the angle ϕ, while the radial coordinate denotes the azimuthal angle θ in the interval [0, π/2] in steps of π/6. c Evolution of the angles θ and ϕ as a function of g_eff for ℏΩ = 10 meV. d Nearest neighbor spin–spin correlation functions $\langle {S}_{i}^{z}{S}_{j}^{z}\rangle$. e Expectation value of the plaquette flux operator W_p. f Average photon number ${n}_{{{{\rm{av}}}}}=\langle \hat{n}\rangle$ of the photo ground state. g Schematic of a fourth-order direct Ru − Ru process illustrating the 1/Ω enhancement of J at low photon energies. h Schematic of a fourth-order ligand-mediated process illustrating the 1/Δ_pd cut-off of K in the low-frequency limit. In all panels, the light-matter coupling refers to the effective single-mode coupling g_eff.

To find how the magnetic ground state flows through the phase diagram of the extended Kitaev model as a function of light-matter coupling and cavity frequency (see Fig. 2b), the magnetic interactions are parameterized by $\bar{J}=\sin \theta \cos \phi$, ${\bar{K}}=\sin \theta \sin \phi$ and $\bar{\Gamma }=\cos \theta$. Here $\bar{J}$, $\bar{K}$, and $\bar{\Gamma }$ are the cavity-renormalized spin parameters divided by the energy $E=\sqrt{{J}^{2}+{K}^{2}+{\Gamma }^{2}}$. The spin parameters for g_eff = 0 were obtained via first principles calculations as discussed in the “Methods” section, and give ϕ = −89.43^∘ and θ = 69.44^∘ placing the equilibrium system in a zigzag antiferromagnetic state adjacent to the ferromagnetic Kitaev QSL. This is in good agreement with experimental findings^43,44. To identify the magnetic state we calculate the spin–spin correlation function $\langle {S}_{i}^{z}{S}_{j}^{z}\rangle$ and the expectation value of the plaquette flux operator ${W}_{p}={2}^{-6}{S}_{i}^{x}{S}_{j}^{y}{S}_{k}^{z}{S}_{l}^{x}{S}_{m}^{y}{S}_{n}^{z}$. In W_p the site indexes traverses a given hexagon p in the positive direction, and the spin component at a given site is determined by the bond pointing out from the hexagon. Magnetic phase transitions are identified by the second-order derivatives of the total energy.

Figure 2b shows the paths traced out by the magnetic ground state as g_eff is increased from the initial value g_eff,i = 0 to the final value g_eff,f = 0.5, for a number of photon energies in the range ℏΩ ∈ [0, 20] meV. We note that parts of the magnetic phase diagram may not be reachable in α-RuCl₃, where the spin parameters are connected via the underlying hopping amplitudes. However, the paths traced out in Figs. 2b and 3a correspond to physical states as they arise from a microscopic electronic Hamiltonian. For all parameters the system flows away from the ferromagnetic Kitaev point, and depending on ℏΩ either remain in the zigzag state or enter the ferromagnetic domain. In particular, for photon energies below ℏΩ ≈ 10 meV, the PGS evolves from a zigzag antiferromagnetic state into a ferromagnetic state as the light-matter coupling is increased (see Fig. 2d and e).

The microscopic origin of the zigzag antiferromagnetic to ferromagnetic phase transition can be traced back to an enhancement of the ratio J/K of the Heisenberg exchange to the Kitaev interaction (see the normalized spin parameters of the zero photon sector displayed in Supplementary Figure 3). In particular, the Kitaev interaction changes sign as a function of light-matter coupling at a point coinciding with the magnetic phase transition, and is thus small in the surrounding region. The increase in J arises from a class of fourth-order electronic processes where a single charge excitation moves between Ru ions while simultaneously emitting or absorbing a virtual cavity photon. One such process is illustrated in Fig. 2g, where the steps (2) → (3) and (3) → (4) both preserve the charge state and thereby lead to an energy difference ~Ω between the virtual states. Since the magnetic interactions are inversely proportional to the energy difference between subsequent virtual states, this leads to an enhancement J ~ 1/Ω in the low-frequency limit. This behavior has been validated in small electron-photon clusters (see Supplementary Figure 2), which show good qualitative agreement with the larger-scale spin-photon model simulations.

The situation is different for the Kitaev interaction, for which the dominant contribution comes from virtual electronic processes proceeding over the ligands. For such processes, there is no way to fourth order in the hoppings to preserve the charge state between two subsequent virtual states (see Fig. 2h). Therefore, the energy difference between participating virtual states has a lower cut-off ~ Δ_pd + Ω. In the low-frequency limit Ω → 0 this leads to a divergence in the ratio J/K ~ (Δ_pd + Ω)/Ω driving the transition to the ferromagnetic state. We note that although the divergence in J will be cut off in any realistic system, such that for Ω < κ (with κ the cavity decay rate) we should have J ~ 1/(Ω + κ), for a cavity with quality factor Q ~ 1000 this still allows for an enhancement of J down to frequencies of at least κ ~ 0.1 meV.

The mechanism driving the magnetic transition is thus consistent with a simultaneous suppression of the effective ligand-mediated hopping and an increase of direct Ru–Ru hopping. However, on physical grounds we also speculate that large-scale fluctuation of the photon field, enhanced in the low-frequency limit, couple indirectly to the ferromagnetic order parameter 〈M〉 = ∑_ij〈S_i ⋅ S_j〉. Such a coupling could potentially be even further enhanced in the macroscopic limit, and could possibly be investigated with quantum Monte Carlo methods⁵¹.

Magnetic phase diagram of a pumped cavity

The results of Fig. 2 show that quantum fluctuations of the cavity vacuum are sufficient to bring about a change in the magnetic ground state of α-RuCl₃ in the THz regime. In contrast, for photon energies large compared to the magnetic energy scale (K/ℏΩ ≪ 1), photon fluctuations are strongly suppressed and the effects of the cavity vacuum are negligible. As will be illustrated below, the situation is radically different when the cavity contains real photons. Due to the suppression of fluctuations the Hamiltonian in the high-frequency regime is well approximated by projecting it onto definite photon sectors, resulting in the diagonal approximation $H\approx {\sum }_{{{{\bf{n}}}}}{{{{\mathcal{H}}}}}_{s,{{{\bf{nn}}}}}\left\vert {{{\bf{n}}}}\right\rangle \left\langle {{{\bf{n}}}}\right\vert$. For large photon energies spin and photon fluctuations therefore decouple, and the magnetic ground state can be obtained by solving the pure spin model $H={{{{\mathcal{H}}}}}_{s,{{{\bf{nn}}}}}$ in each photon sector with effective cavity-renormalized parameters.

The description of the bright cavity in terms of an effective spin model corresponds to a non-equilibrium experimental situation, since the cavity contains real photons, and is therefore referred to as a driven cavity. Assuming that the photon state is quasi-stationary, this approach allows to address the effects of having real photons in the cavity without performing an explicit time evolution. This approximate description of the driven cavity will be valid on short time-scales compared to the inverse cavity decay rate κ⁻¹, and for a cavity quality factor on the order of Q = ℏΩ/κ ~ 1000 or larger the magnetic phase transitions predicted at optical frequencies should be observable. A more detailed discussion of the relation to explicit pumping protocols is provided in the “Methods” section below.

The magnetic phase diagram of α-RuCl₃, interacting with a single effective cavity mode with average photon occupation n_av = 1 and circular polarization, is shown in Fig. 3. The results are displayed as a function of the effective light-matter coupling ${\bar{g}}_{{{{\rm{eff}}}}}=2{g}_{{{{\rm{eff}}}}}\sqrt{{n}_{{{{\rm{av}}}}}}$, obtained by varying g_eff over the interval [0, 0.5] while keeping n_av = 1 fixed, and again the cavity is seen to make the magnetic ground state flow towards the antiferromagnetic Kitaev point (Fig. 3a). The cavity has a drastic effect on the spin parameters as the photon energy is tuned through the charge resonances of the underlying electronic system, as illustrated by the sharp variations of the spin–spin correlation function and plaquette flux around ℏΩ = 2 eV (see Fig. 3c and d as well as Supplementary Figure 4). For a given light-matter coupling this leads to a series of magnetic phase transitions as a function of frequency, stabilizing quasi-stationary states with either ferromagnetic, antiferromagnetic, or spiral order depending on the parameters.

Fig. 3: Magnetic phases in a seeded cavity.

a Paths through the magnetic phase diagram traced out by the spin system as the light-matter coupling ${\bar{g}}_{{{{\rm{eff}}}}}$ increases from ${\bar{g}}_{{{{\rm{eff}}}},i}=0$ to ${\bar{g}}_{{{{\rm{eff}}}},f}=0.5$. The paths are colored according to the cavity photon energy ℏΩ, and the small arrows show the direction of the paths. The polar coordinate denotes the angle ϕ, while the radial coordinate denotes the azimuthal angle θ in the interval [0, π/2] in steps of π/6. b Evolution of the angles θ and ϕ as a function of g for ℏΩ = 3 eV. c Spin–spin correlation function $\langle {S}_{i}^{z}{S}_{j}^{z}\rangle$ as a function of photon energy and light-matter coupling. d Average of the plaquette flux operator W_p as a function of photon energy and light-matter coupling. e Schematic of a third-order process where electrons accumulate a net phase around isosceles triangles, leading to an induced magnetic field. In all panels, the light-matter coupling refers to the effective single-mode coupling ${\bar{g}}_{{{{\rm{eff}}}}}=2{g}_{{{{\rm{eff}}}}}\sqrt{{n}_{{{{\rm{av}}}}}}$ with n_av = 1.

In addition, there is a large region of light-matter couplings ${\bar{g}}_{{{{\rm{eff}}}}}=0.2-0.6$ and photon energies ℏΩ = 2.8 − 3.2 eV where the system enters the antiferromagnetic Kitaev QSL state. The presence of a Kitaev QSL has been verified by calculations of the spin–spin correlation function as well as the plaquette flux operator W_p, which takes on a quantized value W_p = ± 1 in the QSL phase⁵. To ascertain our identification of the QSL phase we also calculated the ground state degeneracy (see Supplementary Figure 5) of the magnetic system, which is forced to be four in the gapped QSL phase in the macroscopic limit⁵². In line with previous work⁵³ we find a robust ground state degeneracy over the whole region identified as the QSL, strongly suggesting that the region identified as the QSL is a gapped and topological quantum spin liquid.

Although the equilibrium system is found to have a ferromagnetic Kitaev interaction, the QSL state appears in a region of antiferromagnetic Kitaev interaction implying a cavity-mediated sign inversion of K. The sign change of K is likely driven by the cavity frequency becoming larger than some charge excitations of the electronic system, which would change the sign of energy denominators of the form ~(E−ℏΩ)⁻¹. In addition, a magnetic field is found to arise from third-order processes involving two Ru–Cl and one Ru–Ru hoppings, leading to a net accumulated flux through the Ru–Cl–Ru triangles (see Fig. 3e). In equilibrium the two fluxes ϕ_R and ϕ_L, arising from right-handed and left-handed hopping paths respectively, are equal and opposite in sign leading to a net cancellation and a resulting zero magnetic field. However, in presence of the chiral cavity field the degeneracy between the right-handed and left-handed hopping paths is broken, and B is in general non-zero.

The QSL is found to be stabilized in a region corresponding to a cavity-induced increase in the magnetic parameters K and B, as well as a simultaneous decrease of J and Γ. A likely explanation for these trends is that in this parameter region the cavity produces a relative enhancement of ligand-mediated hopping as compared to direct Ru–Ru hopping, which would increase K and B relative to the other magnetic parameters. This in consistent with ℏΩ becoming larger than the Ru charge excitations with energies E = 1.8–2.2 eV and approaching the Cl charge-transfer excitations at around 4 eV.

As seen from Fig. 3, α-RuCl₃ also enters the QSL phase in a narrow region of cavity mode frequencies around ℏΩ = 2.1 eV. This coincides with a region where the effective exchange interaction J changes sign, and therefore is very small. More specifically, since the photo-modified magnetic exchange is given close to a resonance by J ~ 1/(E − Ω), and two Ru charge resonances appear close together at about ℏΩ = 2.0 eV and ℏΩ = 2.2 eV, the different sign of J above the lower and below the upper resonance implies a sign reversal between the resonances. This gives a relative amplification of the Kitaev interaction K resulting in a QSL state.

Discussion

The magnetic state of α-RuCl₃ as modified by an optical cavity can be interrogated by various optical techniques. In particular, Kerr and Faraday rotation measurements probe a non-zero magnetization (and therefore the ferromagnetic state)^1,2,9, while linear dichroism measurements are sensitive to the zigzag antiferromagnetic order⁵⁴. While for the QSL state there unfortunately do not exist any direct probes, indirect evidence can be gather via the quasi-particle excitation spectrum as measure by Raman spectroscopy⁵⁵, or measurements of the quantized transverse thermal conductivity⁵⁶.

The results of Figs. 2 and 3 show that the magnetic state of α-RuCl₃ is sensitive to the interaction with a cavity field, and depending on the photon frequency, average cavity occupation and strength of the light-matter coupling, the (quasi-) equilibrium magnetic order can be transformed into any of the magnetic states supported by the extended Kitaev model of Eq. (1). In particular, our results show that the cavity vacuum fluctuations alone are indeed sufficient to change the magnetic order of the photo ground state when the light-matter coupling is sufficiently strong (but within the range of present experimental light-matter couplings). However, as the value of the light-matter coupling is strongly dependent on the cavity size, and scales as $g \sim 1/\sqrt{A}$ with the area of the system, it seems as though the effect of the cavity vanishes quickly in the thermodynamic limit. This naive scaling argument might fail for several reasons, the most important being the neglect of additional cavity modes with momenta q ≈ 0. In fact, the density of modes with momenta smaller than some given cut-off q_m scales as $\rho \sim \sqrt{A}$^24,36,57, and can potentially cancel the volume scaling of g. For the effective single-mode approximation employed here, the mode volume appearing in the light-matter coupling should therefore be interpreted as an effective unit cell mode volume to be fixed either by additional calculations or by comparison to experiment^23,24,45.

An additional subtlety of the thermodynamic limit is the potential coupling of the cavity modes to collective fluctuations of the magnetic system. In particular, close to a magnetic phase boundary, collective fluctuations with a correlation length ξ ~ L (where L is the linear size of the system) are expected to appear. Such collective modes can be expected to play a crucial role in stabilizing the magnetic order, and to couple more strongly to the cavity field⁵¹. However, an explicit treatment of such macroscopic effects, while simultaneously retaining an exact description of spin and photon quantum fluctuations, requires a more advanced methodology going beyond the current work.

For practical purposes, we note that for a given light-matter coupling g_eff the effective coupling ${\bar{g}}_{{{{\rm{eff}}}}}$ can be enhanced by increasing the photon occupation of the cavity. In particular, since ${\bar{g}}_{{{{\rm{eff}}}}}=2g\sqrt{{n}_{{{{\rm{av}}}}}}$, a substantial effective coupling (and corresponding modification of the magnetic interactions in the spin-photon Hamiltonian of Eq. (1)) can be obtained with relatively few photons. Although the present work only includes cavity dissipation on a phenomenological level, the qualitative conclusions are expected to be robust against such effects on time scales short compared to the inverse cavity decay rate κ⁻¹. This includes in particular the identification of a zigzag antiferromagnetic to ferromagnetic transition induced by the cavity vacuum fluctuations, and the presence of a large Kitaev QSL domain in the magnetic phase diagram in the few photon regime. Together these results demonstrate the feasibility of c-QED engineering of the magnetic state of the proximate spin liquid α-RuCl₃ in the few photon limit. Looking ahead, we note that the cavity photon momentum constitutes a tuning parameter in addition to the light-matter coupling, and we therefore envisage that c-QED can be used to also modify phases with finite-q ordering vectors (such as charge density waves). Our work thereby paves the way for utilizing c-QED to control exotic magnetic states in real cavity-embedded quantum materials^20,23.

Methods

Effective low-energy lattice electron Hamiltonian

The crystal structure of a single layer α-RuCl₃ consists of a hexagonal lattice of Ru atoms, surrounded by edge-sharing octahedra of Cl atoms. Due to the strong octahedral crystal field the Ru³⁺d-orbitals are split into a lower e_g and a higher t_2g manifold, the latter consisting of the orbitals d_xy, d_xz, and d_yz. Due to the strong spin-orbit coupling the t_2g states are further split into a lower j_eff = 3/2 quartet, fully occupied in the ground state, and a j_eff = 1/2 doublet with a single hole. The Cl⁻ ions are assumed to be completely filled in the ground state, with their p_x, p_y, and p_z orbitals separated from the d-orbitals by a charge-transfer energy Δ_pd.

The local part of the Hamiltonian is given by^38,39

$$\begin{array}{ll}{H}_{U}\,=\,U\mathop{\sum}\limits_{i\alpha }{\hat{n}}_{i\alpha \uparrow }{\hat{n}}_{i\alpha \downarrow }+\mathop{\sum}\limits_{i\sigma {\sigma }^{{\prime} },\alpha < \beta }({U}^{{\prime} }-{J}_{H}{\delta }_{\sigma {\sigma }^{{\prime} }}){\hat{n}}_{i\alpha \sigma }{\hat{n}}_{i\beta {\sigma }^{{\prime} }}\\ \qquad\quad+\,{J}_{H}\mathop{\sum}\limits_{i,\alpha \ne \beta }({\hat{c}}_{i\alpha \uparrow }^{{\dagger} }{\hat{c}}_{i\alpha \downarrow }^{{\dagger} }{\hat{c}}_{i\beta \downarrow }{\hat{c}}_{i\beta \uparrow }-{\hat{c}}_{i\alpha \uparrow }^{{\dagger} }{\hat{c}}_{i\alpha \downarrow }{\hat{c}}_{i\beta \downarrow }^{{\dagger} }{\hat{c}}_{i\beta \uparrow })\\ \qquad\quad+\,{\Delta }_{pd}\mathop{\sum}\limits_{{i}^{{\prime} }\sigma }{\hat{n}}_{i\sigma }^{p}+\frac{\lambda }{2}\mathop{\sum}\limits_{i}{\hat{{{{\bf{c}}}}}}_{i}^{{\dagger} }({{{\bf{L}}}}\cdot {{{\bf{s}}}}){\hat{{{{\bf{c}}}}}}_{i}.\end{array}$$

(3)

Here U is the intra-orbital Hubbard interaction, ${U}^{{\prime} }$ the inter-orbital interaction, and J_H the Hund’s coupling between the Ru d-orbitals α, β ∈ {yz, xz, xy}. Further, Δ_pd is the single-particle charge transfer energy to add a hole to the Cl ligands, λ is the strength of the spin-orbit coupling (SOC), and the vector of operators is ${\hat{{{{\bf{c}}}}}}_{i}^{{\dagger} }=({\hat{c}}_{iyz\uparrow }^{{\dagger} },{\hat{c}}_{iyz\downarrow }^{{\dagger} },{\hat{c}}_{ixz\uparrow }^{{\dagger} },{\hat{c}}_{ixz\downarrow }^{{\dagger} },{\hat{c}}_{ixy\uparrow }^{{\dagger} },{\hat{c}}_{ixy\downarrow }^{{\dagger} })$. In this basis, the inner product of orbital and spin angular momentum may be written

$${{{\bf{L}}}}\cdot {{{\bf{s}}}}=\left(\begin{array}{ccc}0&-i{\sigma }_{z}&i{\sigma }_{y}\\ i{\sigma }_{z}&0&-i{\sigma }_{x}\\ -i{\sigma }_{y}&i{\sigma }_{x}&0\end{array}\right).$$

(4)

To account for hopping processes between the Ru d-orbitals and Cl p-orbitals, we assume the system has C₃ symmetry. The hopping processes along a Z-bond are described by

$$\begin{array}{ll}{H}_{t}^{{\prime} }\,=\,\mathop{\sum}\limits_{\langle ij\rangle \sigma }\left[\left({\hat{c}}_{iyz\sigma }^{{\dagger} }\,{\hat{c}}_{ixz\sigma }^{{\dagger} }\,{\hat{c}}_{ixy\sigma }^{{\dagger} }\right)\left(\begin{array}{ccc}{t}_{1}&{t}_{2}&{t}_{4}\\ {t}_{2}&{t}_{1}&{t}_{4}\\ {t}_{4}&{t}_{4}&{t}_{3}\end{array}\right)\right.\left(\begin{array}{c}{\hat{c}}_{jyz\sigma }\\ {\hat{c}}_{jxz\sigma }\\ {\hat{c}}_{jxy\sigma }\end{array}\right)\\ \qquad\quad+\,{t}_{pd}\left({\hat{p}}_{{i}^{{\prime} }\sigma }^{{\dagger} }{\hat{c}}_{ixz\sigma }+{\hat{c}}_{jyz\sigma }^{{\dagger} }{\hat{p}}_{{i}^{{\prime} }\sigma }\right.\\ \qquad\quad+\,\left.\left.{\hat{p}}_{{j}^{{\prime} }\sigma }^{{\dagger} }{\hat{c}}_{iyz\sigma }+{\hat{c}}_{jxz\sigma }^{{\dagger} }{\hat{p}}_{{j}^{{\prime} }\sigma }\right)+H.c.\right].\end{array}$$

(5)

Here, t_i with i ∈ {1, 2, 3, 4} parameterize direct d–d hopping processes, while t_pd determines the strength of ligand-mediated hopping via the Cl p_z-orbitals.

Light-matter coupling

In presence of the cavity the total Hamiltonian in the dipole approximation is given by $H={H}_{U}+{H}_{t}+{\sum }_{\lambda }\hslash {\Omega }_{\lambda }{\hat{n}}_{\lambda }$. Here, ℏΩ_λ is the photon energy of a mode described by the label λ and ${\hat{n}}_{\lambda }$ is the corresponding number operator. Inside the cavity the Hamiltonian H_t is modified by the replacements ${\hat{c}}_{i\alpha \sigma }^{{\dagger} }{\hat{c}}_{i\beta \sigma }\to {e}^{i{\phi }_{ij}}{\hat{c}}_{i\alpha \sigma }^{{\dagger} }{\hat{c}}_{i\beta \sigma }$, where the Peierls phases are ${\phi }_{ij}=(ea/\hslash ){{{{\bf{d}}}}}_{ij}\cdot \hat{{{{\bf{A}}}}}$, d_ij = r_j − r_i is the vector between atomic sites i and j measured in units of the Ru–Ru distance a, and the quantum vector potential is

$$\hat{{{{\bf{A}}}}}=\mathop{\sum}\limits_{\lambda }({A}_{\lambda }{{{{\bf{e}}}}}_{\lambda }{\hat{a}}_{\lambda }^{{\dagger} }+{A}_{\lambda }^{* }{{{{\bf{e}}}}}_{\lambda }^{* }{\hat{a}}_{\lambda }).$$

(6)

For a given mode λ this defines the light-matter coupling ${g}_{\lambda }=(ea/\hslash ){A}_{\lambda }=ea/\sqrt{2{\epsilon }_{0}\hslash {\Omega }_{\lambda }V}$. Thus, the Peierls phases can be written as ${\phi }_{ij}={{{{\bf{d}}}}}_{ij}\cdot \hat{{{{\bf{A}}}}}$ with

$$\hat{{{{\bf{A}}}}}=\mathop{\sum}\limits_{\lambda }({g}_{\lambda }{{{{\bf{e}}}}}_{\lambda }{\hat{a}}_{\lambda }^{{\dagger} }+{g}_{\lambda }^{* }{{{{\bf{e}}}}}_{\lambda }^{* }{\hat{a}}_{\lambda }),$$

(7)

and the modified hopping Hamiltonian is

$$\begin{array}{ll}{H}_{t}\,=\,\mathop{\sum}\limits_{\langle ij\rangle \sigma }\left[\left({\hat{c}}_{iyz\sigma }^{{\dagger} }\,{\hat{c}}_{ixz\sigma }^{{\dagger} }\,{\hat{c}}_{ixy\sigma }^{{\dagger} }\right){e}^{i{\phi }_{ij}}\left(\begin{array}{ccc}{t}_{1}&{t}_{2}&{t}_{4}\\ {t}_{2}&{t}_{1}&{t}_{4}\\ {t}_{4}&{t}_{4}&{t}_{3}\end{array}\right)\left(\begin{array}{c}{\hat{c}}_{jyz\sigma }\\ {\hat{c}}_{jxz\sigma }\\ {\hat{c}}_{jxy\sigma }\end{array}\right)\right.\\ \qquad\quad+\,{t}_{pd}\left({e}^{i{\phi }_{{i}^{{\prime} }i}}{\hat{p}}_{{i}^{{\prime} }\sigma }^{{\dagger} }{\hat{c}}_{ixz\sigma }+{e}^{i{\phi }_{j{i}^{{\prime} }}}{\hat{c}}_{jyz\sigma }^{{\dagger} }{\hat{p}}_{{i}^{{\prime} }\sigma }\right.\\ \qquad\quad+\,\left.\left.{e}^{i{\phi }_{{j}^{{\prime} }i}}{\hat{p}}_{{j}^{{\prime} }\sigma }^{{\dagger} }{\hat{c}}_{iyz\sigma }+{e}^{i{\phi }_{j{j}^{{\prime} }}}{\hat{c}}_{jxz\sigma }^{{\dagger} }{\hat{p}}_{{j}^{{\prime} }\sigma }\right)+H.c.\right].\end{array}$$

(8)

Hamiltonian in the photon number basis

The Hamiltonian is expanded in the photon number basis $\left\vert {{{\bf{n}}}}\right\rangle =\left\vert {n}_{1},{n}_{2},\ldots ,{n}_{N}\right\rangle$ according to^30,58,59

$$\begin{array}{ll}H\,=\,\mathop{\sum}\limits_{{{{\bf{nm}}}}}({{{{\bf{1}}}}}_{e}\otimes \left\vert {{{\bf{n}}}}\right\rangle \left\langle {{{\bf{n}}}}\right\vert )H({{{{\bf{1}}}}}_{e}\otimes \left\vert {{{\bf{m}}}}\right\rangle \left\langle {{{\bf{m}}}}\right\vert )\\ \quad\,\,=\,\mathop{\sum}\limits_{{{{\bf{nm}}}}}{H}_{{{{\bf{nm}}}}}\otimes \left\vert {{{\bf{n}}}}\right\rangle \left\langle {{{\bf{m}}}}\right\vert .\end{array}$$

(9)

Here 1_e is the identity operator in the electronic Hilbert space, and the Hamiltonian H_nm is given by

$${H}_{{{{\bf{nm}}}}}=\left({H}_{U}+\mathop{\sum}\limits_{\lambda }\hslash {\Omega }_{\lambda }{n}_{\lambda }\right){\delta }_{{{{\bf{nm}}}}}+{H}_{t,{{{\bf{nm}}}}}.$$

(10)

The matrix elements $\left\langle {{{\bf{n}}}}\right\vert {e}^{i{{{{\bf{d}}}}}_{ij}\cdot \hat{{{{\bf{A}}}}}}\left\vert {{{\bf{m}}}}\right\rangle$ are obtained by noting that since $[{\hat{a}}_{\lambda }^{{\dagger} },{\hat{a}}_{{\lambda }^{{\prime} }}]=0$, the Peierls phases factorize over different modes and

$$\left\langle {{{\bf{n}}}}\right\vert {e}^{i{{{{\bf{d}}}}}_{ij}\cdot \hat{{{{\bf{A}}}}}}\left\vert {{{\bf{m}}}}\right\rangle =\mathop{\prod}\limits_{\lambda }\left\langle {n}_{\lambda }\right\vert {e}^{i{{{{\bf{d}}}}}_{ij}\cdot {\hat{{{{\bf{A}}}}}}_{\lambda }}\left\vert {m}_{\lambda }\right\rangle ,$$

(11)

where ${\hat{{{{\bf{A}}}}}}_{\lambda }={g}_{\lambda }{{{{\bf{e}}}}}_{\lambda }{\hat{a}}_{\lambda }^{{\dagger} }+{g}_{\lambda }^{* }{{{{\bf{e}}}}}_{\lambda }^{* }{\hat{a}}_{\lambda }$.

The single-mode expressions are calculated by introducing the variables η_ijλ = g_λ(d_ij ⋅ e_λ). Using the Baker-Hausdorff formula to expand the exponential, the matrix elements are $\left\langle {n}_{\lambda }\right\vert {e}^{i{{{{\bf{d}}}}}_{ij}\cdot {\hat{{{{\bf{A}}}}}}_{\lambda }}\left\vert {m}_{\lambda }\right\rangle ={i}^{| {n}_{\lambda }-{m}_{\lambda }| }{j}_{{n}_{\lambda },{m}_{\lambda }}^{ij}$ where

$$\begin{array}{lll}{j}_{{n}_{\lambda },{m}_{\lambda }}^{ij}\,=\,{e}^{-| {\eta }_{ij\lambda }{| }^{2}/2}\qquad\qquad\qquad\quad{{\mbox{if}}}\,{n}_{\lambda }\ge {m}_{\lambda }\\ \qquad\qquad\times \,\mathop{\sum }\limits_{k=0}^{{m}_{\lambda }}\frac{{(-1)}^{k}| {\eta }_{ij\lambda }{| }^{2k}{({\eta }_{ij\lambda })}^{{n}_{\lambda }-{m}_{\lambda }}}{k!(k\,+\,{n}_{\lambda }\,-\,{m}_{\lambda })!}\sqrt{\frac{{n}_{\lambda }!}{{m}_{\lambda }!}}\frac{{m}_{\lambda }!}{({m}_{\lambda }\,-\,k)!}\end{array}$$

(12)

$$\begin{array}{ll}{j}_{{n}_{\lambda },{m}_{\lambda }}^{ij}\,=\,{e}^{-| {\eta }_{ij\lambda }{| }^{2}/2}\qquad\qquad\qquad{{\mbox{if}}}\,{m}_{\lambda }\ge {n}_{\lambda }\\ \qquad\qquad\quad\times \,\mathop{\sum }\limits_{k=0}^{{n}_{\lambda }}\frac{{(-1)}^{k}| {\eta }_{ij\lambda }{| }^{2k}{({\eta }_{ij\lambda }^{* })}^{{m}_{\lambda }-{n}_{\lambda }}}{k!(k\,+\,{m}_{\lambda }\,-\,{n}_{\lambda })!}\sqrt{\frac{{m}_{\lambda }!}{{n}_{\lambda }!}}\frac{{n}_{\lambda }!}{({n}_{\lambda }\,-\,k)!}.\end{array}$$

(13)

Using the notation g = {g₁, g₂, …, g_N} to denote the full set of light-matter couplings, the light-matter coupling is described by the function

$${J}_{{{{\bf{nm}}}}}^{ij}({{{\bf{g}}}})=\left\langle {{{\bf{n}}}}\right\vert {e}^{i{{{{\bf{d}}}}}_{ij}\cdot \hat{{{{\bf{A}}}}}}\left\vert {{{\bf{m}}}}\right\rangle =\mathop{\prod}\limits_{\lambda }{j}_{{n}_{\lambda },{m}_{\lambda }}^{ij}.$$

(14)

With this function, the hopping Hamiltonian in presence of the cavity can be written as

$$\begin{array}{ll}{H}_{t,{{{\bf{nm}}}}}\,=\,\mathop{\sum}\limits_{\langle ij\rangle \sigma }\left[\left({\hat{c}}_{iyz\sigma }^{{\dagger} }\,{\hat{c}}_{ixz\sigma }^{{\dagger} }\,{\hat{c}}_{ixy\sigma }^{{\dagger} }\right)\right.\\ \qquad\qquad\times \,{J}_{{{{\bf{nm}}}}}^{ij}({{{\bf{g}}}})\left(\begin{array}{ccc}{t}_{1}&{t}_{2}&{t}_{4}\\ {t}_{2}&{t}_{1}&{t}_{4}\\ {t}_{4}&{t}_{4}&{t}_{3}\end{array}\right)\left(\begin{array}{c}{\hat{c}}_{iyz\sigma }\\ {\hat{c}}_{ixz\sigma }\\ {\hat{c}}_{ixy\sigma }\end{array}\right)\\ \qquad\qquad+\,{t}_{pd}\left({J}_{{{{\bf{nm}}}}}^{{i}^{{\prime} }i}({{{\bf{g}}}}){\hat{p}}_{{i}^{{\prime} }\sigma }^{{\dagger} }{\hat{c}}_{ixz\sigma }+{J}_{{{{\bf{nm}}}}}^{j{i}^{{\prime} }}({{{\bf{g}}}}){\hat{c}}_{jyz\sigma }^{{\dagger} }{\hat{p}}_{{i}^{{\prime} }\sigma }\right.\\ \qquad\qquad+\,\left.\left.{J}_{{{{\bf{nm}}}}}^{{j}^{{\prime} }i}({{{\bf{g}}}}){\hat{p}}_{{j}^{{\prime} }\sigma }^{{\dagger} }{\hat{c}}_{iyz\sigma }+{J}_{{{{\bf{nm}}}}}^{j{j}^{{\prime} }}({{{\bf{g}}}}){\hat{c}}_{jxz\sigma }^{{\dagger} }{\hat{p}}_{{j}^{{\prime} }\sigma }\right)+H.c.\right].\end{array}$$

(15)

Single-mode light-matter coupling and polarization dependence

For a single photon mode the index λ can be dropped, and the function ${J}_{{{{\bf{nm}}}}}^{ij}({{{\bf{g}}}})={J}_{nm}^{ij}(g)={j}_{n,m}^{ij}$. Writing the intersite vector on a given bond as ${{{{\bf{d}}}}}_{ij}=d(\cos {\theta }_{ij}\hat{{{{\bf{x}}}}}+\sin {\theta }_{ij}\hat{{{{\bf{y}}}}})$, the parameter η_ij is given for a circularly polarized field by ${\eta }_{ij}=g({{{{\bf{d}}}}}_{ij}\cdot {{{\bf{e}}}})=gd{e}^{\pm i{\theta }_{ij}}$, where ± denotes a left-handed/right-handed polarization. Noting that ∣η_ij∣ = ∣η∣ = gd so that the spatial dependence of the function ${j}_{n,m}^{ij}$ factorizes, we have ${j}_{n,m}^{ij}={j}_{n,m}{e}^{\pm i(n-m){\theta }_{ij}}$ with

$${j}_{n,m}={e}^{-| \eta {| }^{2}/2}\mathop{\sum }\limits_{k=0}^{m}\frac{| \eta {| }^{2k+n-m}}{k!(k+n-m)!}\sqrt{\frac{n!}{m!}}\frac{m!}{(m-k)!}$$

(16)

for n ≥ m and similarly but with n ↔ m for m ≥ n.

For linear polarization the polarization vector can be written as ${{{\bf{e}}}}=\cos \phi \hat{{{{\bf{x}}}}}+\sin \phi \hat{{{{\bf{y}}}}}$, and we have ${\eta }_{ij}=gd\cos ({\theta }_{ij}-\phi )$.

Connection to Floquet Hamiltonian

As a consistency check, the Hamiltonian H_t,nm for a single circularly polarized photon mode is shown to reduce to the Floquet Hamiltonian of ref. ³⁸ in the semi-classical limit. For a single photon mode, we can drop the mode index, and the function ${J}_{{{{\bf{nm}}}}}^{ij}({{{\bf{g}}}})={J}_{nm}^{ij}(g)={j}_{n,m}^{ij}$. Taking the polarization vector of the photon mode to be ${{{\bf{e}}}}=(\hat{{{{\bf{x}}}}}\pm i\hat{{{{\bf{y}}}}})/\sqrt{2}$, and writing ${{{{\bf{d}}}}}_{ij}=d(\cos {\theta }_{ij}\hat{{{{\bf{x}}}}}+\sin {\theta }_{ij}\hat{{{{\bf{y}}}}})$, we have ${\eta }_{ij}=g({{{{\bf{d}}}}}_{ij}\cdot {{{\bf{e}}}})=gd{e}^{\pm i{\theta }_{ij}}$.

To establish the connection to the Floquet Hamiltonian, what remains is now to show that the function j_n,m(g) reduces to the Bessel function J_n−m(A) in the semi-classical limit. To do this, we write n = m + l and take the limit m → ∞ and g → 0 while keeping $A=2gd\sqrt{m}$ constant. Noting that

$$\begin{array}{ll}{j}_{m+l,m}\,=\,{e}^{-| \eta {| }^{2}/2}\mathop{\sum }\limits_{k=0}^{m}\frac{| \eta {| }^{2k+l}}{k!(k\,+\,l)!}\sqrt{\frac{(m\,+\,l)!}{m!}}\frac{m!}{(m\,-\,k)!}\\ \qquad\quad\,=\,{e}^{-{A}^{2}/8n}\mathop{\sum }\limits_{k=0}^{m}\frac{{\left(\frac{A}{2}\right)}^{2k+l}}{k!(k\,+\,l)!}\sqrt{\frac{(m\,+\,l)!}{m!{m}^{l}}}\frac{m!}{(m\,-\,k)!{m}^{k}},\end{array}$$

(17)

and that (m + l)!/m! → m^l and m/(m − k)! → m^k as m → ∞, we have

$$\mathop{\lim }\limits_{m\to \infty }{j}_{m+l,m}=\mathop{\sum }\limits_{k=0}^{\infty }\frac{{(-1)}^{k}{\left(\frac{A}{2}\right)}^{2k+l}}{k!(k+l)!}={J}_{l}(A).$$

(18)

Since l = n − m the Hamiltonian H_t,nm reproduces the correct Floquet Hamiltonian in the semi-classical limit.

Cavity-modified spin Hamiltonian

The total Hamiltonian H is down-folded to the spin sector by using quasi-degenerate perturbation theory to eliminate H_t order by order⁶⁰. This perturbation expansion is performed numerically up to fourth order in t/U using exact diagonalization of the Ru₂Cl₂ cluster shown in Fig. 3e. This allows to extract the nearest-neighbor magnetic interactions and their dependence on the light-matter coupling, cavity frequency and polarization. The structure of the perturbation expansion allows the down-folding to be performed separately within each photon sector, giving a coupled spin-photon Hamiltonian of the form

$${{{\mathcal{H}}}}=\mathop{\sum}\limits_{{{{\bf{nm}}}}}\left({{{{\mathcal{H}}}}}_{s,{{{\bf{nm}}}}}+{\delta }_{{{{\bf{nm}}}}}\mathop{\sum}\limits_{\lambda }\hslash {\Omega }_{\lambda }{n}_{\lambda }\right)\left\vert {{{\bf{n}}}}\right\rangle \left\langle {{{\bf{m}}}}\right\vert .$$

(19)

Here ${{{{\mathcal{H}}}}}_{s,{{{\bf{nm}}}}}$ is the spin Hamiltonian in the photon sector connecting photon numbers n and m, and is given by

$$\begin{array}{ll}{{{{\mathcal{H}}}}}_{s,{{{\bf{nm}}}}}\,=\,\mathop{\sum}\limits_{\langle ij\rangle }\left(\begin{array}{ccc}{S}_{i}^{\alpha }&{S}_{i}^{\beta }&{S}_{i}^{\gamma }\end{array}\right){\left(\begin{array}{ccc}J&\Gamma &{\Gamma }^{{\prime} }\\ \Gamma &J&{\Gamma }^{{\prime} }\\ {\Gamma }^{{\prime} }&{\Gamma }^{{\prime} }&J+K\end{array}\right)}_{{{{\bf{nm}}}}}\left(\begin{array}{c}{S}_{j}^{\alpha }\\ {S}_{j}^{\beta }\\ {S}_{j}^{\gamma }\end{array}\right)\\ \qquad\qquad+\,{B}_{{{{\bf{nm}}}}}\mathop{\sum}\limits_{i}{\hat{{{{\bf{e}}}}}}_{B}\cdot {{{{\bf{S}}}}}_{i}\end{array}$$

(20)

where each bond (ij) is labeled by the indexes αβ(γ) ∈ {xy(z), yz(x), zx(y)}. In this Hamiltonian, the magnetic parameters J, K, and Γ all depend on the light-matter coupling g and the photon numbers n and m.

Exact diagonalization

The ground state of the coupled spin-photon system was obtained by exact diagonalization of a 24-site spin cluster interacting with a single photon mode. To perform the calculations we used to open-source Python package QuSpin^61,62.

Number and coherent states in the driven cavity

In the discussions above the driven cavity is approximately described by projecting the full spin-photon Hamiltonian onto a definite photon number sector. Here, we provide additional details of the relation between this approach and an explicit cavity driving protocol.

The cavity can be driven to a finite photon population through a laser-cavity interaction ${H}_{{{{\rm{pump}}}}}=f(t)\sin (\Omega t)({a}^{{\dagger} }+a)$, where f is some envelope function. Assuming that the driving is performed in a time T short compared to the characteristic magnetic time scale K⁻¹, the photons will at time T be in a coherent state with some average occupation n_av. Since the cavity contains real photons, the coupled spin-photon problem now corresponds to a non-equilibrium situation, and should be explicitly treated using open quantum systems methods. Since this is computationally unfeasible if an exact description of the correlated spins is to be retained, we assume that the photon state is quasi-stationary (which will hold over some time-scale κ⁻¹ set by the cavity decay rate) and study the effects of having real photons in the cavity in this quasi-equilibrium setting.

The difference between having photons in a coherent state or in a number state is expected to have a negligible influence on the magnetic properties, and we therefore choose the later approach for computational convenience. This can be understood by recalling that in the high-frequency limit we neglect the coupling between different photon number sectors (i.e., there is no cross-coupling), which implies that the expectation value calculated in a coherent state can be obtained by weighting the expectation values 〈A〉_n obtained in different photon number sectors by their projection on the coherent state. More precisely, we have that the expectation value of an operator A evaluated in the coherent state is related to the expectation values in the number states by

$$\langle A\rangle =\mathop{\sum}\limits_{n}{e}^{-| \alpha {| }^{2}}\frac{{\alpha }^{2n}}{n!}{\langle A\rangle }_{n}.$$

(21)

Since 〈A〉_n has a much weaker dependence on n than the prefactor, we can write $\langle A\rangle \approx (\,{{\mbox{const}}}\,)\times {\langle A\rangle }_{{n}_{{{{\rm{av}}}}}}$ showing that these expectation values are proportional.

Model parameters from first principles

To determine the parameters of the electronic Hamiltonians H_U and H_t, we performed first principles simulations of monolayer α-RuCl₃ with the Octopus electronic structure code^63,64. The single particle parameters where obtained from a Wannierization of the Ru 4d and Cl 3p orbitals in the paramagnetic state using Wannier90^65,66, while the interaction parameters where determined using the hybrid DFT+U functional ACBN0 in the zigzag state⁶³. We have checked that the single particle parameters differ by less than one percent between the paramagnetic and zigzag states. The resulting electronic parameters are given in Table 1, the equilibrium spin parameters in Table 2, and the electronic band structure is shown in Fig. 4.

Table 1 Electronic parameters of α-RuCl₃.

Table 2 Equilibrium spin parameters of α-RuCl₃.

Fig. 4: Electronic band structure of α-RuCl₃.

Highest valence bands (red) and lowest conductions bands (blue) of α-RuCl₃ obtained within the DFT+U formalism⁶³.

The calculations were performed in a $1\times \sqrt{3}$ supercell to account for the zigzag magnetic structure, using the experimental lattice parameters a = 5.98 Å and b = 10.35 Å. Mixed boundary conditions, periodic in the in-plane direction and open in the out-of-plane direction, where used together with a vacuum region of 15 Å to ensure convergence in the out-of-plane direction. A 8 × 8k-point grid and a real-space grid spacing of 0.3 Bohr were employed. Using the ACBN0 functional a self-consistent effective interaction U_eff = U − J_H was determined on both the Ru and Cl ions, and the Kanamori parameters U and J_H where calculated in the final states after convergence had been reached.

For a two-dimensional cavity the lowest energy photon mode has a frequency Ω = πc/L_z, where L_z is the extent of the cavity in the z-direction. Therefore, the light-matter coupling is independent of the photon frequency, and is given by

$$g=\frac{ea}{\sqrt{2{\epsilon }_{0}\hslash \Omega {a}_{x}{a}_{y}{L}_{z}}}=\frac{ea}{\sqrt{2\pi {\epsilon }_{0}\hslash c{a}_{x}{a}_{y}}}=0.12$$

(22)

assuming a_xa_y = a².