Engineering and Probing Non-Abelian Chiral Spin Liquids Using Periodically Driven Ultracold Atoms

We propose a scheme to implement Kitaev’s honeycomb model with cold atoms, based on a periodic (Floquet) drive, in view of realizing and probing non-Abelian chiral spin liquids using quantum simulators. We derive the effective Hamiltonian to leading order in the inverse-frequency expansion, and show that the drive opens up a topological gap in the spectrum without mixing the effective Majorana and vortex degrees of freedom. We address the challenge of probing the physics of Majorana fermions, while having access only to the original composite spin degrees of freedom. Specifically, we propose to detect the properties of the chiral spin liquid phase using gap spectroscopy and edge quenches in the presence of the Floquet drive. The resulting chiral edge signal, which relates to the thermal Hall effect associated with neutral Majorana currents, is found to be robust for realistically prepared states. By combining strong interactions with Floquet engineering, our work paves the way for future studies of non-Abelian excitations and quantized thermal transport using quantum simulators.

Popular Summary

Chiral spin liquids are one of the most intriguing quantum phases of matter. Predicted to exist in two-dimensional quantum systems, these topological phases exhibit excitations known as non-Abelian anyons, which go beyond the conventional boson and fermion paradigm. Despite intense efforts in condensed-matter physics, it is still debated whether quantum spin liquid order actually exists in nature. One emblematic instance is the so-called Kitaev honeycomb model, which admits an exact analytical solution. Remarkably, recent advances in the design of quantum simulators open a possible path for the first experimental realization of the original Kitaev honeycomb model, hence suggesting that chiral spin liquids (including their exotic excitations) can be studied and manipulated in a highly controlled experimental environment.

In this work, we propose a realization of the Kitaev honeycomb model using quantum simulators based on nonequilibrium periodic drives. The model exhibits a chiral spin liquid ground state, which features non-Abelian excitations consisting of vortex (flux) and Majorana fractional excitations. We introduce practical methods for probing the topological properties of this exotic state of matter, using well-designed drive sequences and available measurement schemes. In particular, our detection protocols allow us to unambiguously reveal the topological heat current associated with Majorana edge modes, a hallmark signature of chiral spin liquids.

This work paves the way for the quantum simulation of chiral spin liquids, offering an appealing alternative to their experimental investigation in quantum materials.

Figure 1

Model and drive. (a) The Kitaev honeycomb model is an emblematic instance of a ${\mathbb{Z}}_2$ lattice gauge theory; it can be mapped exactly to Majorana fermions defined on the vertices of the honeycomb lattice (cyan) and dispersionless vortices residing on the plaquettes (orange). (b) Dispersion relation of the Majorana fermions (cyan bands) with gap closing at the Dirac cones. In our driving protocol, the vortex dispersion (orange) remains flat, which renders vortices immobile. (c) Ground state phase diagram of the original Kitaev model, featuring three disconnected gapped phases $A_{\alpha }$ (Abelian topological phases equivalent to the toric code) and a gapless spin liquid phase $B$. (d) Floquet realization: periodic modulation of the $x$, $y$, and $z$ bonds at frequency ${\omega }_D$ breaks time-reversal symmetry and opens up a topological gap $\mathrm{\Delta }$ in the Majorana dispersion in phase $B$, leaving the vortex dispersion intact; see (b). The ground state of the effective Hamiltonian $H_{\mathrm{eff}}$ is a chiral spin liquid, a topological phase exhibiting non-Abelian excitations.

Figure 2

Gap spectroscopy. (a) Majorana spectrum, $E_{\mathrm{Majorana}}$, showing a bulk gap opening at the Dirac points as well as chiral edge modes for a system with periodic boundary conditions in the $y$ direction (cylindrical geometry). (b) Many-body energy spectrum, $E_{\mathrm{spin}}$, of the spin Hamiltonian $H_{\mathrm{eff}}$ in the vortex-free sector for a cylindrical geometry (inset). The many-body ground state $|G⟩$ of the system is a product state of the filled lower Majorana band and the ground state of the vortex sector. The bulk many-body gap $\mathrm{\Delta }$ is indicated, and the edge states are located within the corresponding shaded area. (c) The number of excited particles, $N_{\mathrm{exc}}(t_{\ast })$, shows well-defined resonances as a function of the probe frequency ${\omega }_p$ (vertical dotted line), from which the gap can be extracted: $\mathrm{\Delta }={\omega }_{\mathrm{res}}/2$. The inset shows $N_{\mathrm{exc}}$ as a function of time for ${\omega }_p/J=0.23$; $t_{\ast }$ denotes the measurement time ($J{t}_{\ast }=300$). Spectroscopy is performed in the presence of the Floquet drive on an ${18}_x\times 6_y$ cylinder for $A_p=0.02$ and ${\omega }_D/J=12$. (d) The extracted gap decreases when the drive frequency ${\omega }_D$ increases, and shows good agreement with the theoretical value of the bulk gap, $\mathrm{\Delta }=\pi \sqrt{3}{J}^2/4{\omega }_D$, of the static model $H_{\mathrm{eff}}$ (black line). Circles (crosses) correspond to a cylindrical (planar) geometry. The legend shows the number of unit cells used in the calculations.

Figure 3

Detecting chiral Majorana edge transport. (a) Schematic representation of the nonequilibrium quench protocol, which we use to reveal the chiral heat current in the ground state of $H_{\mathrm{eff}}$. (b) The Fourier spectrum associated with $zz$ correlations on the edge, $C_m(\ell )$ [Eq. (7)], which we use to probe the chiral dispersion relation of Majorana excitations. The slope in $(\omega ,k)$ space (dashed line) equals half the group velocity of the edge-mode excitations ($v_{\mathrm{edge}}$). Upper (cyan) and lower (red-hot) panels correspond to positive and negative circularly polarized Floquet drives: $A^{\circlearrowright }(k,\omega )$ and $A^{\circlearrowleft }(k,\omega )$. The system size is ${16}_x\times {32}_y$. The model parameters are ${\omega }_D/J=10$, $J_q/J=0.008$, and $N_T=1001$. (c) The group velocity of the chiral mode $v_{\mathrm{edge}}$ extracted from (b) decreases as one increases the Floquet drive frequency ${\omega }_D$, in agreement with Floquet theory (which predicts a vanishing edge-mode velocity in the infinite-frequency limit). The evolution times are proportional to the drive period, $JT{N}_T\approx 100\times 2\pi$ for every ${\omega }_D$ point. The error bars show the fit error; see Appendix pp9. (d) Same as (b) but for a smaller system of size $8_x\times 8_y$; the remaining parameters are the same as (b).

Figure 4

Detecting chiral Majorana edge transport after quasiadiabatic quantum state preparation. (a) The bulk Majorana gap at the Dirac point closes at the critical point between phases $A_z$ and $B$ during the ramp $J_{\alpha }(t)=Jt/t_{\mathrm{ramp}}$, $\alpha =x,y$ [inset, see also Fig. 1]. The gray area marks the existence of edge modes in the chiral spin liquid phase. (b) Same as Fig. 3 but with the initial state for the quench $U_F\to {\tilde{U}}_F$ prepared from the ground state at $J_x=J_y=0$ with use of the ramp protocol in (a). The gray region at $k=0$ shows a strong signal, disregarded in the postprocessing procedure. Chiral velocity fits (straight dashed lines) are performed in the regions $(k>0,\omega >0{)}_{\circlearrowright }$ and $(k>0,\omega <0{)}_{\circlearrowleft }$. The parameters are $t_{\mathrm{ramp}}=250J$, $N_T=1001$, and ${\omega }_D=10J$. (c) The chiral current velocity $v_{\mathrm{edge}}$, extracted from the fits in (b), as a function of the observation cycle number $N_T$ after the quench (details can be found in Appendix pp9). Longer observation times allow us to better resolve the signal and reduce the fit error. The parameters are $J{t}_{\mathrm{ramp}}=250$ and ${\omega }_D=8J$. (d) Chiral mode velocity fits $v_{\mathrm{edge}}$ as a function of the ramp duration $t_{\mathrm{ramp}}$ for two drive frequencies, ${\omega }_D/J=8$ and ${\omega }_D/J=15$. Horizontal dashed lines correspond to the theoretical values. The number of observation cycles $N_T=$ 801 for ${\omega }_D/J=8$ and $N_T=$ 1501 for ${\omega }_D/J=15$. The system size is $8_x\times {16}_y$ for all panels.

Figure 5

Physical implementation. (a) Short $\pi /2$ pulses in between the steps of the Floquet drive, applied on a timescale much shorter than the inverse coupling strength $1/J$, turn the $z$ interactions into $x$ and $y$ interactions. (b) Fermion implementation of the Floquet drive: (i) Fermi-Hubbard model with periodic modulation of the hopping; see Eq. (12). In the high-frequency limit, $J_{\uparrow }\ll U\ll {\omega }_D$, application of the inverse-frequency expansion (i)$\to$(ii) leads to a modified static fermion Hamiltonian with nearest-neighbor hopping terms (also in the spin channel) of strength $J_{\uparrow }/3$. A Schrieffer-Wolff transformation (iii) then gives the effective spin-$1/2$ Hamiltonian $H_{\mathrm{eff}}$ from Eq. (4) with interaction strength $J=2J_{\uparrow }^2/9U$, which governs the physics of the lower Hubbard band stroboscopically. Alternatively, (i)$\to$(iv), in the strongly interacting limit $J_{\uparrow }\ll {\omega }_D\ll U$, we first apply the Schrieffer-Wolff transformation to arrive at the spin-$1/2$ Hamiltonian (13), with periodic-in-time, spatially oscillating $z$ interaction $J_{zz,ij}(t)=2J_{\uparrow ,ij}^2(t)/U$. A subsequent application of the inverse-frequency expansion then gives rise to $H_{\mathrm{eff}}$ with interaction strength $J=2J_{\uparrow }^2/3U$. Step (iv) can also be independently taken as the starting point for the implementation in a spin-1/2 quantum simulator (red-shaded oval). (c) Time evolution of the $z$ and $x$ magnetization on a single honeycomb cell (open boundary conditions) for a few different values of the drive frequency ${\omega }_D/J_{\uparrow }$. The initial states are product states along the spin $z$ and $x$ directions (see the legend). The fermionic Floquet implementation (solid lines) agrees well with the dynamics generated by the spin-only Hamiltonian $H_{\mathrm{eff}}$ (dotted lines) over a wide range of drive frequencies ${\omega }_D$, except in the vicinity of resonances where $l{\omega }_D\approx U$ ($l\in \mathbb{N}$). The model parameters are $J_{\downarrow }/J_{\uparrow }=0$ and $U/J_{\uparrow }=20$, with $J_z(\zeta )$ given in Eq. (14).

Figure 6

Schematic representation of the notation $[ijk{]}_{\alpha \beta \gamma }$, showing the three-body terms to order ${\omega }_D^{-1}$, which are responsible for opening the topological bulk gap. The letters refer to the spin operator applied to each vertex, and the bond colors correspond to the bonds in the Kitaev honeycomb model; see Fig 1.

Figure 7

Numerical test of the effective Hamiltonian and kick operator against an exact simulation of a spin-$1/2$ system. Here $U_{\mathrm{eff}}^{(0)}=\exp (-iT{H}_{\mathrm{eff}}^{(0)})$ and $U_{\mathrm{eff}}^{(0+1)}=\exp (-i{K}_{\mathrm{eff}}^{(1)}(0))\exp (-iT[H_{\mathrm{eff}}^{(0)}+H_{\mathrm{eff}}^{(1)}])\exp (i{K}_{\mathrm{eff}}^{(1)}(0))$. The data show clearly that all corrections to the desired order are captured by Eq. (4). Note that the transformation generated by the kick operator is crucial for the validity of the test. The system size is a single honeycomb cell, and all interactions are isotropic, $J_{\alpha }^{\prime }=J^{\prime }$.

Figure 8

Schematic representation of the notation used in Eq. (B2). The $p\mathrm{th}$ plaquette $p$ contains the six lattice sites $p_1$–$p_6$, with $b_1$ and $b_2$ denoting the two $z$ bonds on the plaquette. The unit cell contains sites from two sublattices marked by open and filled circles. $x$, $y$, and $z$ refer to the direction of the bounds in the Hamiltonian.

Figure 9

The number of excited particles depends on the probe frequency ${\omega }_p$ for both planar geometry and cylindrical geometry with different system sizes (see the legend). (a) For the planar geometry with a small system size of $6_x\times 6_y$, the boundary states have a non-negligible influence, which shows up as an excitation of a midgap state. (b) For a large system size, ${18}_x\times {18}_y$, the influence of the boundary states becomes negligible, and the lowest resonance peak indicates the location of the bulk gap. (c),(d) For a cylindrical geometry, the lowest resonance peak indicates the location of the bulk gap quite well for the ${12}_x\times 6_y$ system and the ${18}_x\times 6_y$ system. The vertical dotted black line indicates the lowest resonant frequency ${\omega }_{\mathrm{res}}$, and the dashed green line indicates the location of $2\mathrm{\Delta }$. The insets show the time dependence of the signal $N_{\mathrm{exc}}(t)$ at ${\omega }_p={\omega }_{\mathrm{res}}$ used to extract the bulk gap. The other parameters are the same as for Fig. 2, i.e., $J{t}_{\ast }=300$, $A_p=0.02$, and ${\omega }_D/J=12$.

Figure 10

Schematic representation of the notation for the quench protocol: $m$ labels a unit cell at the edge of the system; the site $i$ belongs to the unit cell $m$, while $⟨i,j_x{⟩}_x$, $⟨i,j_y{⟩}_y$, and $⟨i,j_z{⟩}_z$ indicate the bonds connected to $i$ that are considered in the expression for the local unit cell energy $H_m$; see Eq. (6).

Figure 11

Fourier spectrum $|A(k,\omega )|$ of the local unit cell energy and the $zz$ correlations following the boundary quench in a ${16}_x\times {32}_y$ system. (a),(b) Fourier spectrum for the local unit cell energy for ${\omega }_D/J=8$ and $N_T=801$ and for ${\omega }_D/J=20$ and $N_T=2001$, respectively. The chiral signal is clearly noticeable among additional excitations. (c),(d) Fourier spectrum for the $zz$ correlations for ${\omega }_D/J=8$ and $N_T=801$ and for ${\omega }_D/J=50$, and $N_T=5001$, respectively. The dominant chiral behavior is clearly discernible, and the amount of nonchiral excitations appears reduced. For all panels, the slope of the dashed black line corresponds to the theoretical prediction, $2v_{\mathrm{edge}}^{\mathrm{th}}\times k+0.5$, and $J_q/J=0.008$. For comparison, in Fig. 3, the dashed black line corresponds to a fit (see Appendix pp9).

Figure 12

Fourier spectrum $|A(k,\omega )|$ of the $zz$ correlations at the boundary of a ${16}_x\times {32}_y$ system following the quench (see the main text) for two different points belonging to phase $A_z$ of the phase diagram [Fig. 1]. (a) At the critical point ($J_x=J_y=0.5J$ and $J_z=J$) the chiral bulk gap in the Majorana spectrum closes [see Fig. 4], and the chiral signal is lost. (b) The chiral signal is not observable deeper in the toric code phase ($J_x=J_y=0.25J$ and $J_z=J$) either. The figure demonstrates that the chiral signal in the Fourier spectrum is an intrinsic property of the chiral spin liquid phase $B$; in particular, it is not induced by the “polarization” of the Floquet drive. The parameters are ${\omega }_D/J=8$, $J_q/J=0.008´$, and $N_T=801$, and the remaining parameters are the same as for Fig. 11.

Figure 13

The effect of increasing quench strength $J_q$ on the Fourier spectrum of the $zz$ correlation and corresponding edge-state velocity. The plots show the Fourier spectrum of the $zz$ correlation for three different quench strengths: (a) $J_q/J=0.008$, (b) $J_q/J=0.16$, and (c) $J_q/J=0.6$. The corresponding extracted edge-state velocity $v_{\mathrm{edge}}$ is shown in (d). As the quench strength increases, residual, “noisy” components appear in the Fourier spectrum, but the dominant chiral signal remains qualitatively unchanged. The system size is ${16}_x\times {32}_y$ for a cylindrical geometry. The model parameters used in the study are ${\omega }_D/J=10$ and $N_T=1001$, and the cutoff used in the fitting is 0.4.

Figure 14

Fourier spectrum of the $zz$ correlation for various ramp and observation times. (a) For a short ramp time, $J{t}_{\mathrm{ramp}}=250$, a strong signal appears at $k=0$, which dominates the chiral signal shown in Fig. 4. (b) At long ramp times, $J{t}_{\mathrm{ramp}}=5000$, the dominance of the chiral signal is restored. For (a),(b), $N_T=1001$ produces a good resolution of the signal. (c) For short observation times, $N_T=101$, the resolution becomes low but the chiral behavior still exists. The drive frequency ${\omega }_D/J=10$ for (a)–(c). (d) For large drive frequencies (e.g., ${\omega }_D/J=50$), a strong signal at $k=0$ comparable to the chiral signal persists even at long ramp times $J{t}_{\mathrm{ramp}}=5000$. The system size is $8_x\times {16}_y$ for all panels, and the slope of the dashed black line corresponds to the theoretical prediction: $2v_{\mathrm{edge}}^{\mathrm{th}}$.

Figure 15

Fourier spectrum of the $zz$ correlation calculated by (a) the linear response method and (b) the exact numerical simulation. The system size is ${16}_x\times {32}_y$ and the model parameters are ${\omega }_D=8/J$, $J_q/J=0.008$, and $N_T=800$. The dashed black line corresponds to the theoretical prediction: $2v_{\mathrm{edge}}^{\mathrm{th}}k+1/2$; see Eq. (H16).

Figure 16

Edge-state velocity $v_{\mathrm{edge}}$ extracted from fits to $|A(k,\omega )|$ with different cutoff values at $J{t}_{\mathrm{ramp}}=250$. For a small cutoff value of $0.2$, residual excitations dominate the Fourier signal, which leads to larger error bars. Additionally, for larger values of ${\omega }_D/J$, the gap is smaller, and the strength of the chiral signal is comparable to the strength of residual excitations due to the finite ramp and observation times, which also results in larger error bars even if the cutoff value is relatively large. Here $N_T=50{\omega }_D/J+1$ is chosen long enough to resolve the chiral signal. The system size is $8_x\times {16}_y$ (cylindrical geometry).