Emergent strong zero mode through local Floquet engineering

Periodically driven quantum systems host exotic phenomena that often do not have any analog in undriven systems. Floquet prethermalization and dynamical freezing of certain observables, via the emergence of conservation laws, are realized by controlling the drive frequency. Recent experimental developments in synthetic quantum matter, such as superconducting qubits and cold atoms, have opened avenues for implementing local Floquet engineering which can achieve spatially modulated quantum control of states. Here, we uncover the novel memory effects of local periodic driving in a nonintegrable spin-half staggered Heisenberg chain. For a boundary-driven protocol at the dynamical freezing frequency, we show the formation of an approximate strong zero mode, a prethermal quasilocal operator, due to the emergence of a discrete global ${\mathbb{Z}}_2$ symmetry. This is captured by constructing an accurate effective Floquet Hamiltonian using a higher-order partially resummed Floquet-Magnus expansion. The lifetime of the boundary spin can be exponentially enhanced by enlarging the set of suitably chosen driven sites. We demonstrate that in the asymptotic limit, achieved by increasing the number of driven sites, a strong zero mode emerges, where the lifetime of the boundary spin grows exponentially with system size. The nonlocal processes in the Floquet Hamiltonian play a pivotal role in the total freezing of the boundary spin in the thermodynamic limit. The novel dynamics of the boundary spin is accompanied by a rich structure of entanglement in the Floquet eigenstates where specific bipartitions yield an area-law scaling while the entanglement for random bipartitions scales as a volume-law. Our work addresses the long-standing question of the existence of a strong zero mode in a nonintegrable model and elucidates the complex nature of thermalization in locally driven systems.

Figure 1

A schematic figure showing the model interaction and drive protocol adapted in this work. Red (blue) color is used to denote ferro- (antiferro-) magnetic exchange in consecutive links. A stationary global field is applied in the $z$ direction, whereas a time-periodic field in the $x$ direction is applied at a few specific sites.

Figure 2

Integrated level spacing distribution $I\left(s\right)$ for $H_0$ (in blue) and $H_0$ without staggering (in yellow), both with $h=1,N=16$. Black lines show a comparison with Poisson statistics, the signature of integrability (dashed), and GOE, valid for chaotic systems (dotted).

Figure 3

Stroboscopic dynamics of (a) single site entanglement entropy $S_i^{\text{vN}}$ for all the sites ($i$) for $N=14$. The red (dashed) line is used to denote the maximum single-site entropy ($ln2$). $\overline{S_1^{\text{vN}}}$ (averaged $S_1^{\text{vN}}$ over first 2000 cycles) is plotted with $\omega$ in the inset showing the absence of dynamic freezing of the edge spin even for the most appropriate choice of the drive frequency. A vertical dashed line is also drawn at ${\omega }_c^1=7.5$ to clearly show the minima in $\overline{S_1^{\text{vN}}}$ (${\omega }_m^1$) is shifted from ${\omega }_c^1$. (b) entanglement entropy of the boundary site ($S_1^{\text{vN}}$) for different system sizes ($N$). $h=1$, $\gamma =15$, and $\omega ={\omega }_m^1=7.53$ for both the panels. The $N=16$ and $N=18$ calculations were done with matrix-product-state-based TEBD calculations for a time step of $\delta t=\tau /100$ and a maximum bound dimension of $2^{N/2}$.

Figure 4

Comparison of (a) real, (b) imaginary, and (c) absolute values of $H_F(2,1)$ between exact numerics (black dots) and expressions obtained from F-M expansions. Blue (green) curve is obtained from Eq. (18) [Eq. (19)]. A vertical red (dashed) line is used in (c) to guide the eye to the minima of $H_F(2,1)$ which is clearly shifted from ${\omega }_c^1=7.5$. $h=1$, $\gamma =15$, and $N=8$.

Figure 5

Comparison of stroboscopic dynamics from exact numerics and analytical $H_F$ obtained from the resummed F-M expansion at (a) $\omega =5$ (b) $\omega =7.53$ (c) $\omega =10$. $\gamma =15$, $h=1$, and $N=8$. Away from the freezing frequency the stroboscopic dynamics is well captured for all spins within a “local” approximation for $H_F$, in comparison to exact numerics. At or near the freezing frequency, nonlocal terms are required to accurately capture the behavior of the driven boundary spin.

Figure 6

Stroboscopic $S_i^{\text{vN}}$ (top) and $\langle S_i^x\rangle$ (bottom) for different two-site driving. Locations of the driven sites ($s_d$) are mentioned at the top of each column. $h=1,\gamma =15,\omega ={\omega }_m^1=7.53,\text{and}N=14$.

Figure 7

Long time stroboscopic dynamics of $S_1^{\text{vN}}$ for different multi-site driving showing the optimal choice of drive locations (every fourth site) to freeze the boundary spin. Inset shows $S_1^{\text{vN}}$ vs $n$ for $s_d=\{1,5\}$ for different system sizes. All parameters are the same as in Fig. 6.

Figure 8

Emergence of strong zero-mode: (a) $\overline{|\langle D_z\rangle |}$ vs $\omega$ showing the emergence of a discrete symmetry, $D_z={\prod }_i\left(2S_i^z\right)$, at $\omega \sim {\omega }_m^1$. (b) $\overline{|\langle {\mathrm{\Phi }}^{+}|S_1^x|{\mathrm{\Phi }}^{-}\rangle |}$ vs $\omega$. Inset shows that the condition $S_1^x|{\mathrm{\Phi }}^{\pm }\rangle =\frac{1}{2}|{\mathrm{\Phi }}^{\mp }\rangle$ is satisfied better for $s_d=\{1,5\}$ and $s_d=\{1,5,9\}$ compared to $s_d=\left\{1\right\}$. (c) $\overline{\left|\mathrm{\Delta }\right|}$ vs $\omega$ showing a deep at $\omega \sim {\omega }_m^1$ which implies the emergence of SZM induced quasidegeneracy between the Floquet eigenvalues from different $D_z$ sectors. $N=10$, $\gamma =15$ for all the plots. We have performed a moving average over the raw data (in an $\omega$ interval of 0.008) for panel (b). The vertical dashed line is a guide to the eye to the position of strongest freezing.

Figure 9

Scaling of averaged pairing gaps ($\overline{\left|\mathrm{\Delta }\right|}$) (a) with system size ($N$) for a fixed number of driven sites and (b) with increasing the size of the driven cluster $n_d^{\text{4th}}$ for a fixed system size ($N=14$). All other parameters are the same as in Fig. 6. The dashed green line in subfigure (b) is a best-fit curve to an exponential ansatz, $\overline{\left|\mathrm{\Delta }\right|}\sim exp(-2.9n_d^{\text{4th}})$.

Figure 10

Entanglement entropy of the boundary site with the rest of the system ($S_1^{\text{vN}}$) in the upper panels, and sites {1,5,9,13} with the rest of the system ($S_{1,5,9,13}^{\text{vN}}$) in the lower panels for $N=14$. The driven sites are mentioned at the top of each column. All other parameters are the same as in Fig. 6. Page values of entanglement entropies are shown by the red dashed line. The bipartition used to calculate the entanglement entropy in (c) and (d) is shown in (e).

Figure 11

The matrix element $H_F(2,1)$ vs $\omega$ for global driving comparing $H_F^{\left(2\right)}$ with exact numerical data. $\gamma =15$, $h=1$, and $N=8$.

Figure 12

Level statistics of the Floquet eigenvalues, (a) at a generic drive frequency $\omega =5.6$ (b) at the freezing frequency $\omega ={\omega }_m^1=7.53$. The blue line corresponds to the Wigner surmise (an approximation to the GOE statistics of chaotic systems), while the red line is the Poisson distribution (integrable systems). $N=14$, $\gamma =15$.

Figure 13

The behavior of $|\langle {\mathrm{\Phi }}_r^{-}|S_1^x|{\mathrm{\Phi }}_r^{+}\rangle |$ in the upper panels (a-d) in parallel with the behaviors of ${\mathcal{R}}_r$ in the lower panels (e)–(h). The driven sites are mentioned on top of the upper panels (valid for both panels in each column). The values of $\overline{|\langle {\mathrm{\Phi }}^{-}|S_1^x|{\mathrm{\Phi }}^{+}\rangle |}$ are (a) 0.4275, (b) 0.4942, (c) 0.4996, (d) 0.4998 and the values of $\overline{\mathcal{R}}$ are (e) 0.8550, (f) 0.9883, (g) 0.9991, (h) 0.9996. $\gamma =15$, $\omega =7.53$, $N=14$ for all the plots.

Figure 14

The behavior of $S_1^{\text{vN}}$ in all the Floquet eigenstates in the upper panels in parallel with the behavior of $|{\delta }^{\pm }|$ (not distinguished) in the lower panels. The average deviations of $S_1^{\text{vN}}$ from the maximal value $ln2$ (shown in red dashed line) are (a) $1.21\times {10}^{-2}$, (b) $1.42\times {10}^{-4}$, (c) $8.89\times {10}^{-5}$, (d) $8.34\times {10}^{-5}$, and the average values of $|{\delta }^{\pm }|$ are (e) 0.11, (f) $1.38\times {10}^{-2}$, (g) $1.09\times {10}^{-2}$, and (h) $1.05\times {10}^{-2}$. All other parameters are the same as in Fig. 13.

Figure 15

Comparison of the scaling of $\overline{\left|\mathrm{\Delta }\right|}$ and $\overline{|{\mathrm{\Delta }}^{\prime }|}$. The latter is calculated from exact numerics as well as using a semianalytical expression [Eq. (F3)]. Parameters: $h=1,\gamma =15,\omega ={\omega }_m^1=7.53,\text{and}N=14$.