Quantum Information Scrambling in Quantum Many-body Scarred Systems

Quantum many-body scarred systems host special non-thermal eigenstates that support periodic revival dynamics and weakly break the ergodicity. Here, we study the quantum information scrambling dynamics in quantum many-body scarred systems, with a focus on the"PXP"model. We use the out-of-time-ordered correlator (OTOC) and Holevo information as measures of the information scrambling, and apply an efficient numerical method based on matrix product operators to compute them up to 41 spins. We find that both the OTOC and Holevo information exhibit a linear light cone and periodic oscillations inside the light cone for initial states within the scarred subspace, which is in sharp contrast to thermal or many-body localized systems. The periodic revivals of OTOCs and Holevo information signify unusual breakdown of quantum chaos and are not equivalent to the revival dynamics of state fidelity or local observables studied in the previous literature. To explain the formation of the linear light cone structure, we provide a perturbation-type calculation based on a phenomenological model. In addition, we demonstrate that the OTOC and Holevo information dynamics of the"PXP"model can be measured using the Rydberg-atom quantum simulators with current experimental technologies, and numerically identify the measurable signatures using experimental parameters.


I. INTRODUCTION
Isolated quantum many-body systems would eventually thermalize under the time evolution and their subsystems relax to the equilibrium, leading to the emergence of ergodicity and statistical mechanics. During this process, any local information preserved in the initial states scrambles into the entire system and becomes unrecoverable. This kind of quantum thermalization phenomena has been illustrated by the eigenstate thermalization hypothesis (ETH) in the past decades [1,2]. While numerous works have confirmed the universality and correctness of ETH in various scenarios [3][4][5][6], discovering quantum many-body systems violating the ETH is still of fundamental and practical importance. Known exceptions to the ETH paradigm include the integrable [7] and manybody localized (MBL) [8,9] systems, which have either exact or approximate extensive conserved quantities to prevent the systems from thermalization. Recently, in experiments with Rydberg atoms, non-thermal periodic revival dynamics has been observed after quenching the system from a special high-energy Néel state [10,11]. Subsequent theoretical works attribute this weak ergodicity breaking to a small fraction of ETH-violating eigenstates (dubbed "quantum manybody scars") embedded in a sea of thermal eigenstates [12][13][14] . Here, we investigate the scrambling of quantum information in many-body scarred systems (see Fig. 1 for a pictorial illustration).
Quantum information scrambling describes the propagation and effective loss of initial local information in quantum many-body dynamics. It has attracted considerable attention in different contexts, including the black hole thermaldynamics [15][16][17], quantum many-body quench dynamics [18][19][20][21], Individual neutral atoms are trapped with optical tweezers (vertical red beams). The global Rydberg laser (the horizontal blue beam), with the Rabi frequency Ω and detuning ∆, couples the atomic ground state |g (black circles) with the Rydberg excited state |r (white circles). A pair of Rydberg atoms at the distance Rij shares the van der Waals repulsion Uij = U0/R 6 ij . (b) A schematic illustration of the information scrambling dynamics in systems of different thermalization classes. and machine learning [22][23][24]. Except for a few pronounced integrable examples [25][26][27][28][29][30][31][32][33][34][35], for general quantum many-body systems with short-range interactions obeying the ETH, the Lieb-Robinson bound [36] restricts the correlation propagation within a linear light cone (analogous to the causal light cone in relativistic theories). In contrast, strong disorders in MBL systems prevent the local transport, resulting in a logarithmic light cone for information spreading [37][38][39][40][41][42] [see Fig.  1(b) for a sketch]. Quantum many-body scars, as a new thermalization class in many-body dynamics, possess the potential to exhibit different information spreading behaviours from the former two cases. Despite previous extensive studies of quantum many-body scars from various perspectives , the exploration of quantum information scrambling in quan-tum many-body scarred systems is still lacking hitherto.
In this paper, we apply the out-of-time-ordered correlator (OTOC) and Holevo information as measures to study this problem. We find that both the OTOC and Holevo information exhibit a linear light cone and periodic revival dynamics inside the light cone for initial states within the scarred subspace, displaying distinct features from the ETH and MBL cases [ Fig. 1(b)]. Moreover, the persistent oscillations will disappear and the information propagation speed will increase once we choose generic high-energy initial states. The periodic oscillations of OTOCs and Holevo information signify unusual breakdown of quantum chaos and persistent backflow of quantum information. Due to the action of interleaved operators, the dynamics of OTOCs and Holevo information actually involve both the scarred subspace and the thermal eigenstate bath. We emphasize that their dynamics can not be readily deduced from the eigenstate decomposition of initial states, thus distinguishing our work from the previous literature [12][13][14]. In order to explain the linear light cone structure, we further provide a perturbation-type calculation in the interaction picture based on a phenomenological model proposed in [43]. Finally, we propose an experiment with Rydberg atoms to observe the predicted exotic OTOC and Holevo information dynamics, and numerically identify the measurable signatures using experimental parameters.

II. MODEL
Motivated by the Rydberg-atom experiment [10], in the limit of Rydberg blockade [72] the physics of quantum manybody scars are extracted as the one-dimensional (1D) "PXP" model with an open boundary condition [12,13]: where P j = (1 − σ z j )/2, σ x,y,z j are Pauli matrices of the j-th qubit, L is the number of total qubits, and | ↓ (↑) represents the atomic ground (Rydberg excited) state |g(r) . Below we always consider the dynamics within the constrained Hilbert space (where computational bases with two nearby up spins | · · · ↑↑ · · · are removed). The PXP Hamiltonian is non-integrable and chaotic according to the level statistics studies, yet it holds a small fraction of ETH-violating scarred eigenstates that support the periodic revival dynamics for initial states within the scarred subspace (such as the Néel state |Z 2 = | ↑↓↑ · · · ↑↓ ). Whereas, for generic high-energy initial states (such as |0 = | ↓↓ · · · ↓ ) the dynamics will quickly become chaotic and no revival occurs [12,13,44].
The OTOC utilizes the Heisenberg operator growth to characterize the information scrambling and quantum chaos, and is defined as [17,73,74] where |ψ is an initial pure state, W i , V j are local observables defined on sites i, j, and V j (t) = e iHt V j e −iHt ( = 1). The OTOC directly connects to the squared commuta- |ψ by the relation C ij (t) = 2(1 − Re(F ij (t))), for unitary operators W i , V j . A simple physical picture for the OTOC is that if the Heisenberg operator growth of V j (t) does not reach the site i, [W i , V j (t)] = 0, F ij (t) = 1, while the equalities will break down when sites i, j become correlated inside the causal region. Note that the OTOCs F ij (t) consist of the forward and backward Hamiltonian evolution, between which there exist interleaved local operators W i and V j . Since the evolved quantum states (e −iHt |ψ or e −iHt W i |ψ ) in general are not the eigenstates of inserting operators W i and V j , the OTOC dynamics can be essentially different from those of simple local observables. For instance, in the study of MBL systems, people have realized the distinction between the dynamics of OTOCs and common two-point correlators: OTOCs exhibit a logarithmic light cone and decay inside [37][38][39][40][41][42], whereas the two-point correlators are bounded to be exponentially small due to the localization nature [8,9]. In the following discussions, we mainly focus on the ZZ-OTOC (W = σ z , V = σ z ) and XZ-OTOC (W = σ x , V = σ z ). The Holevo information (or Holevo χ quantity) originates from the quantum information theory to upper bound the accessible information between two separate agents [75,76]. Consider that Alice prepares mixed states ρ X in the set {ρ 1 , ρ 2 , · · · , ρ n } with probability {p 1 , p 2 , · · · , p n } respectively, and then sends ρ X to Bob. With any kind of positive operator-valued measures (POVMs), the amount of information Bob can obtain about the variable X according to the measurement outcome Y is bounded by is the mutual information between X and Y , S(ρ) = −Tr(ρ log ρ) denotes the von Neumann entanglement entropy. Roughly speaking, the Holevo information describes the distinguishability of states in the set {ρ 1 , ρ 2 , · · · , ρ n }. For instance, p 1 = p 2 = 1/2, if ρ 1 = | ↑ ↑ |, ρ 2 = | ↓ ↓ |, then χ = 1; while if ρ 1 = ρ 2 = | ↑ ↑ |, then χ = 0. Compared with OTOCs, Holevo information possesses richer informationtheoretic meanings and its experimental measurement does not require the inverse Hamiltonian evolution. We hence expect that Holevo information can be used to characterize the information spreading dynamics and solve informationtheoretic problems in more quantum many-body systems [77][78][79].
Here, we regard the reduced Hamiltonian evolution on subsystems as quantum communication channels in the original setup of Holevo information [75] and use it to study the information scrambling dynamics. Consider the Hamiltonian evolution on two different initial states, |ψ and σ x i |ψ , taking |ψ to be a computational basis state (e.g., |Z 2 or |0 ). In other words, initially we have encoded one bit information at the site i [| ↑ (↓) i ]. We demonstrate how the one-bit information scrambles into the entire system by computing the Holevo information on the site j: where ρ j (t) and ρ j (t) are reduced density matrices of the j-th spin after the Hamiltonian evolution for the initial state |ψ and σ x i |ψ respectively. χ j (t) measures how much information one could obtain by any local probe on the j-th site for these two sets of evolution.

III. NUMERICAL SIMULATIONS
In Fig. 2, we numerically compute the OTOC and Holevo information for the PXP model as diagnoses of the information scrambling dynamics. Specifically, we apply an efficient matrix-product-operator (MPO) method [30] to calculate the ZZ-OTOC and XZ-OTOC up to L = 41 spins (see Appendix. C for algorithm details). We observe the following features: For the initial state |ψ = |Z 2 [ Fig. 2(a), (c)], both the ZZ-OTOC and XZ-OTOC spread ballistically, forming a linear light-cone structure with the butterfly velocity v b ∼ 0.6 (inverse of the light cone slope) [16,[80][81][82][83][84]; inside the light cone, the OTOC dynamics show evident periodic revivals with the period T ≈ 4.71, consistent with the oscillation period of state fidelity and local observables in Hamiltonian evolution [12,44]. Besides, the oscillations for different sites j are synchronized, meaning that F ij (t) of different j have the same period T and reach the maxima at the same time t. In contrast, for a generic high-energy initial state like |0 ( Fig. 2(b), (d)), the OTOCs have a larger butterfly velocity v b ∼ 1 (information scrambles faster) and decay without discernible periodic revivals inside the light cone.
We emphasize that except for the ZZ-OTOC of the |Z 2 initial state, the periodic oscillations of OTOCs inside the light cone can not be deduced from the approximate recurrence of the |Z 2 state under Hamiltonian evolution. For instance, in the XZ-OTOC case, the inserting σ x i operators will flip the |Z 2 state partially out of the scarred subspace and introduce the influence from the thermal eigenstate bath. In order to confirm the generality of the revival behaviours, we numerically show that the persistent and synchronized oscillations still appear even if the initial states are replaced with scarred eigenstates. We use the overlap between eigenstates |E n of the PXP Hamiltonian and |Z 2 [ Fig. 2(h)], as well as the halfchain entanglement entropy of |E n (see Appendix. C for details), to distinguish quantum many-body scars from typical thermal eigenstates [12,13]. In Fig. 2(g), we transform the ZZ-OTOC dynamics into the frequency domain and find that: For initial states being scarred eigenstates, there exist ω 0 = 2π/T, 2ω 0 , 3ω 0 peaks in the spectra for different sites j, but not for the case of initial states being the generic thermal eigenstates. Fig. 2(g) indicates that the periodic revivals of OTOCs are general phenomena for initial states within the scarred subspace. The recurrence of quantum information signified by the OTOCs is not equivalent to the recurrence of quantum states through Hamiltonian evolution, for which the energy eigenstates actually have no dynamics. In Appendix. C we further display the OTOC dynamics with initial states being eigenstates of the PXP Hamiltonian and superposition states of |Z 2 and |Z 2 = ( L i=1 σ x i )|Z 2 . Similar information scrambling dynamics also emerge when probed by the Holevo information. As shown in Fig. 2(e), for two sets of Hamiltonian evolution on |Z 2 and σ x L/2 |Z 2 , χ j (t) initially vanishes everywhere except on the central qubit, where one bit local information is encoded. As time evolves, non-zero Holevo information can be probed at other sites j, which forms a linear light cone in spacetime and also persistent oscillations inside the light cone with a period T ≈ 4.73, consistent with the OTOC results. In comparison with the initial states |0 and σ x L/2 |0 [ Fig. 2(f)], χ j (t) peaks when the information wavefront arrives and eventually diminishes for long time t, following the indistinguishability of ρ j (t) and ρ j (t) predicted by the ETH. The information spreading velocity v h for the |Z 2 case is v h ∼ 0.6, less than that of the |0 case v h ∼ 0.9. Similarly, due to the action of the σ x L/2 operator, the Holevo information dynamics involve both the scarred subspace and the thermal eigenstate bath, which are further discussed in Appendix. B. Fig. 2(e) is obtained by the time-evolving block decimation (TEBD) algorithm [85,86] based on the matrix-product-state (MPS) ansatz, which leverages the relatively low entanglement entropy of scarred eigenstates and is not applicable to the |0 case [12,13].

IV. ANALYTICAL EXPLANATIONS
In this section, we provide analytical explanations for the OTOC and Holevo information dynamics observed in numerical simulations. First, we mention that the analytical computation of OTOCs is a challenging task for general quantum many-body systems, despite a few pronounced solvable examples such as the Sachdev-Ye-Kitaev (SYK) model [87][88][89][90]. For the PXP model, it is difficult to directly deal with its OTOC dynamics by the perturbation method since all the terms in the Hamiltonian have the same interaction strength. However, the essential features of the PXP model consist of two parts: the periodic oscillations within and the chaotic dynamics out of the scarred subspace, which we adopt a phenomenological model proposed in [43] to effectively describe: plays the role of thermalization for states out of the scarred subspace. The L + 1 scarred eigenstates are all the x-direction Once we start the Hamiltonian evolution from |ψ = | ↑↑ · · · ↑ , perfect quantum state revivals with period T = 2π/Ω will be observed, imitating most but not all the characteristics of the PXP model [43] (see Appendix. A for model details and relevant numerical results).
For the non-trivial situation of XZ-OTOC W 1 = σ x 1 , V r = σ z r with the initial state |ψ = | ↑↑ · · · ↑ , specially at time points t = nT /2 (n = 0, 1, 2 · · · ), the OTOCs are simplified into F (r, t) = (−1) n φ(t)|σ z r |φ(t) , where |φ(t) = e −iH t σ x 1 |ψ . The physical picture for the forma- tion of linear light cone structure in OTOC and Holevo information dynamics is that: Without the action of σ x 1 , the state |φ(t) will undergo perfect periodic oscillations, thus F (r, t = nT /2) ≡ 1. However, now σ x 1 penetrates the scarred subspace on the first site. The "heat flow" (like a quasi-particle created by the local quench [21]) leaks and propagates ballistically through the entire system, leading to the decay of OTOC F (r, t = nT /2) < 1 after the wavefront of heat flow arrives. Note that in the setup of Holevo information dynamics, we also have a Hamiltonian evolution term e −iHt σ x L/2 |Z 2 , which has the same structure of |φ(t) . We hence deduce that the dynamics of Holevo information and OTOCs are closely related and follow the same physical picture. Indeed the numerical results of the two criteria support each other and exhibit similar information scrambling behaviors.
We use the interaction picture of H 0 to remove the Rabi oscillation effect. For the early growth region of OTOCs, we can split the evolution time t into r pieces with ∆t = t/r, J∆t 1, where J is the average energy scale of all the J µν ij . After some lengthy calculations shown in Appendix. A, the deterioration of perfect oscillations at the site r will be dominantly caused by an operator prod- Here, a is some model-dependent O(1) constant and we have generalized the special time points t = nT /2 to arbitrary t before the wavefront arrives. Eq. (4) depicts a linear light cone t ∝ r/J in the early growth region of OTOCs.
Despite the fact that the OTOC dynamics have linear light cone structure both in the ETH and many-body scarred systems, the thermalization processes are pretty different: In the ETH systems, when we start the Hamiltonian evolution, thermalization happens everywhere globally on the quantum state; On the contrary for many-body scarred systems, according to the calculations above, the periodic revivals are destroyed locally at some sites and the deterioration effect then propagates through the entire system. The global versus local thermalization processes are reminiscent of the spin correlation dynamics in global [20] versus local [21] quench dynamics, and can also explain the larger butterfly velocity for the |0 initial state than that of the |Z 2 initial state in OTOC and Holevo information dynamics.
The perturbation-type calculations break down for the OTOC and Holevo information dynamics deep inside the light cone. As mentioned in the previous section, the persistent and synchronized oscillations inside the light cone still appear even if the initial states are replaced with scarred eigenstates (see Appendix. C Fig. 7 and 9). We attribute the unusual revival dynamics inside the light cone to the local rather than global thermalization argued above and robustness of scarred eigenstates against local perturbations: The quantum many-body scars have been shown to present certain robustness against perturbations and disorders [63][64][65][66][67]. Under some local perturbations like σ x(z) i(j) operators on the evolved quantum states (e −iH t σ x 1 |ψ ), thermalization happens locally and propagates outwards. After the wavefront of thermalization effect passes, the periodic oscillations are preserved to some extent due to the robustness of scar eigenstates, leading to the revival pattern observed in the OTOC and Holevo information dynamics. This picture are also observed in Appendix. A Fig.  4 and Appendix. B Fig. 5, and could be an interesting avenue for future investigations.

V. EXPERIMENT PROPOSAL
Recent experimental progress has enabled the measurements of information scrambling dynamics in well-controlled synthetic quantum systems, including nuclear magnetic resonance systems [91,92], trapped ions [93][94][95][96] and superconducting qubits [97][98][99][100]. The PXP model can be naturally realized with the Rydberg-atom platform [72,101,102], governed by the following Hamiltonian: σ x i connects the atomic ground (|g ) and Rydberg excited (|r ) state on the site i with the Rabi frequency Ω and detuning ∆. Two Rydberg atoms at distance R ij share the van der Waals repulsion approximately reduces to the PXP model, which provides us an opportunity to experimentally measure the ZZ-OTOC and Holevo information dynamics [see Fig. 1(a)].
The main difficulty in simulating OTOC dynamics lies in the implementation of the inverse Hamiltonian evolution exp (−i(−H)t). Fortunately, the PXP model has a particlehole symmetry operation ( i σ z i ) H ( i σ z i ) = −H to reverse the sign of H. For the Rydberg-atom Hamiltonian Eq. (5), the i σ z i operator only changes the sign of the Rabi oscillation term while keeping the detuning term and Rydberg blockade structure intact. Below we denote the Hamiltonians for experimental evolution and inverse evolution as H ± = ±(Ω/2) i σ x i − ∆ i n i + U 1 i n i n i+1 + U 2 i n i n i+2 , where we only take the next-nearest-neighbor interaction U 2 into consideration of error analysis due to the sixth power decay. Specifically, for the case of ZZ-OTOC F ij (t) = ψ|σ z i σ z j (t)σ z i σ z j (t)|ψ and initial states |ψ = |Z 2 or |0 , since σ z i |ψ = (−1) ni+1 |ψ , we have where |Ψ j (t) = e −iH−t σ z j e −iH+t |ψ , and e −iH−t = ( i σ z i ) e −iH+t ( i σ z i ). The measurements of ZZ-OTOC are hence reduced to the evaluation for the expectation value of σ z i on a time-dependent quantum state |Ψ j (t) . The implementation of |Ψ j (t) only requires Hamiltonian evolution of Eq. (5) together with global and individual Pauli-Z gates [72,101]. The measurement protocol for the Holevo information is straightforward: preparing two initial states |Z 2 (0) and σ x L/2 |Z 2 (0) , evolving the system with H independently, and finally doing quantum state tomography for all the single-qubit density matrices to obtain χ j (t) via Eq. (3).
When approximating the PXP model using the lowest U (1) symmetry sector of the Rydberg-atom Hamiltonian, the errors mainly come from the invalidity of the condition U 1 Ω U 2 . Constrained by the geometry of 1D equally-spaced atoms, U 1 /U 2 ∼ 64 is fixed, such that Ω around √ U 1 U 2 can best fulfill the inequality above. Besides, we numerically observe that a small non-zero ∆ in H ± can further eliminate the effect of U 2 , which is probably due to the cancellation of −∆ i n i and U 2 i n i n i+2 terms on the |Z 2 initial state. The fact that non-zero detunings amplify the signatures of quantum many-body scars is reminiscent of the recent experiment [11] in which periodically driven detunings stabilize the scar revivals, and might be of independent research interest. In Fig.  3, we display the numerical simulations of the ZZ-OTOC and Holevo information dynamics for the Rydberg-atom Hamiltonian approximated by H ± , with the initial state |Z 2 , under experimental parameters from [10] and an optimized ∆ to increase the oscillation contrast as much as possible (see Appendix. C for more details). The linear light cone contour and periodic oscillations inside the light cone can be readily observed.

VI. CONCLUSIONS AND DISCUSSIONS
In summary, we have studied the information scrambling dynamics in quantum many-body scarred systems. We found an unconventional paradigm (a linear light cone with periodic oscillations inside, revealed by OTOCs and Holevo information) intrinsically distinct from the previously studied thermal or MBL systems. Based on perturbation-type calculations in the interaction picture, we provided analytical explanations for this paradigm. In addition, we have also proposed an experiment to measure the predicted exotic information scrambling dynamics with current Rydberg atom technologies.
Because of the forward and backward Hamiltonian evolution and interleaved local operators, the periodic revivals of OTOCs and Holevo information inside the light cone can not be directly deduced from the eigenstate decomposition of initial states, thus distinguishing our work from the previous literature. For initial states within the scarred subspace, locally encoded information is retrievable elsewhere even at late time, indicating persistent backflow of quantum infor-mation and unusual breakdown of quantum chaos. Our results present an information-theoretic perspective on quantum many-body scars, which would connect to a number of possible directions for future studies, such as Hilbert space fragmentation [59,60,[103][104][105], the classical OTOC and chaos theory [44-46, 106, 107], robustness of scarred eigenstates under perturbations [63][64][65][66][67], black hole physics [77][78][79], and quantum technology applications including quantum memory and quantum sensing [14,69,108]. In this Appendix, we provide detailed analytical derivations of the OTOC dynamics in the early growth region, i.e. Eq. (4) in the main text, and relevant discussions. As mentioned in the main text, the essential features of the PXP model can be described by a phenomenological model proposed in [43]: These scarred eigenstates form an exact su(2) algebra (while the scars in the PXP model form an approximate one) [43], leading to the periodic revival dynamics governed by the global Rabi oscillation term H 0 , from initial states within the scarred subspace (for example, the z-direction Dicke states like |ψ = | ↑↑ · · · ↑ ). In contrast, initial states out of the scarred subspace can not be projected out by P i,i+1 , thus are affected by the R i,i+3 terms and have chaotic dynamics. The can be viewed as the overlap between two time-dependent quantum states F ij (t) = ψ 2 (t)|ψ 1 (t) : First of all, for the ZZ-OTOC case of the PXP model with the initial state |Z 2 (or the H model above with the initial state |ψ = | ↑↑ · · · ↑ ), since |Z 2 (|ψ ) is the eigenstate of the σ z i(j) operator, the periodic revival dynamics of the quantum state directly give us the periodic and synchronized oscillations of the ZZ-OTOC F ij (t). Note that the oscillation period of the ZZ-OTOC is supposed to be T /2 = π/Ω, which is indeed true for the H model with perfect quantum many-body scars. Because after T /2 evolution, |Z 2 evolves to |Z 2 = ( L i=1 σ x i )|Z 2 (| ↑↑ · · · ↑ evolves to | ↓↓ · · · ↓ ), which is again the eigenstate of the σ z i(j) . However, the periodic revival dynamics and emergent su(2) algebra of the PXP model are not perfect [43], resulting in smaller peak values of F ij (t = (2n + 1)T /2) than those of F ij (t = nT ) (n = 0, 1, 2 · · · ). Eventually we observe an overall oscillation period T = 2π/Ω in numerical simulations.
Specially for the PXP model, since we consider the dynamics within the constrained Hilbert space (where computational bases with two nearby up spins | · · · ↑↑ · · · are removed), only the ZZ-OTOC dynamics for generic initial states |ψ are legal, namely within the constrained Hilbert space. The XZ-OTOCs for the |Z 2 initial state with W i = σ x i acting on an up spin | ↑ i are also legal and have well-defined physical meanings in the overlap interpretation Eq. (A4). Fortunately, the problem of constrained Hilbert space does not appear in the H model above. Besides, the Holevo information dynamics for general initial states are also not well defined, since it will be difficult to locally encode one bit information on general entangled states, for example the scarred eigenstates. Now we consider the non-trivial XZ-OTOC dynamics of H with the initial state |ψ = | ↑↑ · · · ↑ . For W 1 = σ x 1 , V r = σ z r , t = nT /2 (n = 0, 1, 2 · · · ), according to Eq. (A4), Here by utilizing the property of state recurrence, we successfully convert the four-body OTOC into an observable average on a time-dependent quantum state |φ(t) = e −iH t σ x 1 |ψ . In order to further explore the role of the R thermalization term, we adopt the interaction picture of H 0 = (Ω/2) i σ x i to remove the Rabi oscillation effect. We denote the quantum states and operators in the interaction picture with hats: The quantum state evolution in the H 0 interaction picture is where we have introduced the time ordering operator T t . According to Eq. (A6), at the special time points t = nT /2, the OTOCs in the interaction picture have the expression Below we consider the so-called early growth region of the OTOC, namely the time interval from t = 0 to the time t just before the wavefront of heat flow arrives at the site r. Then we can split the evolution time t into r pieces with ∆t = t/r, J∆t 1 (J is the average energy scale of all the J µν ij ), to treat the dynamics perturbatively: whereR According to Eq. (A9) and (A10), besides the leading "1" term in Eq. (A11), the OTOC F (r, t = nT /2) will be dominantly influenced by the following operator product series: The operator product series Eq. (A13) vividly characterizes propagation of the thermalization effect from the site 1 to site r. Since H 0 only contains single-body operators, other R i (n∆t) operator product series acting on σ x 1 |ψ will either be zero due to the projectors P i,i+1 , or unable to reach the site r and affect the perfect revival dynamics there. Inserting Eq. (A13) into Eq. (A10), we obtain that the leading correction term of the OTOC F (r, t = nT /2) can be bounded by where a is some model-dependent O(1) constant, [·] ± denotes commutator(−) and anti-commutator(+) (depending on the parity of r to choose which), and || · || denotes the operator norm. The bound in Eq. (A14) does not depend on T or Ω (which stands for the resolution in time), so we are able to generalize the special time points t = nT /2 to arbitrary t before the wavefront arrives. Finally we have the following OTOC behaviors in the early growth region: Eq. (A15) readily depicts a linear light cone structure t ∝ r/J in the early growth region. While we specifically compute the XZ-OTOC of the phenomenological model Eq. (A1), the physical pictures of local thermalization and ballistic propagation hold for other quantum many-body scarred systems, other OTOCs and Holevo information dynamics, and can be applied to explain the linear light cone structure observed in numerical simulations. An additional remark is that focuses on rotation dynamics for the three middle spins of |Z2 (left) and σ x L/2 |Z2 (right). Due to the kinetic constraints, the retarded spin rotations lead to the periodic distinguishabilty of single-qubit reduced density matrices ρj(t) and ρ j (t). The Holevo information χj(t) hence forms the persistent oscillation pattern inside the light cone. the perturbation calculations above are similar to the derivations of Kubo formula in the linear response theory, and the correction term Eq. (A13) corresponds to the r-th order response to the perturbationR(t) [109].
Numerical simulations for the OTOC dynamics of the H model are displayed in Fig. 4. The persistent and synchronized oscillations are readily observed for initial states within the scarred subspace in plots (a), (b) and (c). In contrast for the Z 2 Néel state, OTOCs constantly decay to zero without any revival. These numerical results further confirm the generality of our conclusions to other models with quantum manybody scars. In future studies it is also interesting to explore whether the formation mechanisms of quantum many-body scars [14] will affect the information scrambling dynamics in corresponding models.
As mentioned in the main text, the persistent and synchronized oscillations inside the light cone still appear even if the initial states are replaced with scarred eigenstates or superposition states of |Z 2 and |Z 2 = ( L i=1 σ x i )|Z 2 (see Fig. 7  and 9), which confirms that periodic revivals of OTOCs are general phenomena for initial states within the scarred subspace, not some fine-tuned results. We have attributed the un-usual revival dynamics inside the light cone to the local rather than global thermalization and robustness of scarred eigenstates against local perturbations. In this sense, the periodic revival behaviors of OTOCs in quantum many-body scarred systems are somewhat similar to the dynamics in time crystals [105,[110][111][112]: The ZZ-OTOC dynamics for the |Z 2 initial state is an analogue to the time crystals without perturbations; while the XZ-OTOCs for |Z 2 and ZZ-OTOCs for general initial states within the scarred subspace correspond to the perturbed time crystals, which will exhibit certain robustness and maintain the synchronized oscillations. This topic might need deeper understanding and more powerful techniques to deal with, so that is expected to inspire more analytical studies in the future. One possible direction is about the relation with the classical OTOC and chaos theory [44-46, 106, 107]. According to the quantum-classical correspondence principle, we may replace the quantum commutators in the OTOC with classical Poisson brackets. By utilizing the semi-classical Lagrangian of quantum many-body scarred systems, like the one proposed in [44] with the timedependent variational principle (TDVP), we are able to compute the Poisson brackets ({·} P.B. ) of some observables like FIG. 6. The time-splitting MPO algorithm for OTOC calculations. The blue squares denote MPS tensors and the white (red) squares denote MPO tensors. All the rectangles stand for the Trotterized Hamiltonian evolution blocks. The total evolution time t has been equally split to each part of Fij(t) = ψ(t/2)|W † i (−t/2)V † j (t/2)Wi(−t/2)Vj(t/2)|ψ(t/2) to reduce the required bond dimension. The overall contraction of the tensor network above will give the OTOC Fij(t).

Appendix B: More about Holevo Information Dynamics
As mentioned in the main text, we regard the reduced Hamiltonian evolution on subsystems as quantum communication channels in the original setup of Holevo information [75] and use it to study the information scrambling dynamics. Through two sets of Hamiltonian evolution on different initial states, |ψ and σ x i |ψ , we demonstrate how onebit local information at the site i scrambles into the entire system by computing the Holevo information on the site j χ j (t) = S (ρ j (t) + ρ j (t))/2 − S(ρ j (t)) + S(ρ j (t)) /2, where ρ j (t) and ρ j (t) are reduced density matrices of the j-th spin after the Hamiltonian evolution for the initial state |ψ and σ x i |ψ respectively. χ j (t) measures how much information one could obtain by any local probe on the j-th site for these two sets of evolution. If the many-body dynamics follow the predictions of ETH and the computational bases |ψ and σ x i |ψ bear the same energy with respect to the Hamiltonian, at the late time, local single-body density matrices ρ j (t) and ρ j (t) will become indistinguishable. The Holevo information χ j (t) will quickly peak when the effect of site i arrives and diminish close to zero afterwards (so does the case of |0 and σ x L/2 |0 ). In contrast, we have observed persistent and synchronized oscillations inside the light cone for the case of |Z 2 and σ x L/2 |Z 2 . When unpacking the dynamics of Holevo information, we find that the kinetic constraints of the PXP model play an important role for the periodic oscillations inside the light cone: According to the PXP Hamiltonian, each spin will rotate freely around the x axis if both its neighbors are in the | ↓ state; otherwise its rotation dynamics are frozen. In Fig.  5, we display all the single-qubit reduced density matrices used for the computation of Holevo information. The density matrix of a spin-1/2 can be written as ρ = (I + a · σ)/2, where σ = (σ x , σ y , σ z ), and a is the Bloch vector on the Bloch sphere, such that the dynamics of ρ(t) can be fully parametrized by a(t). For the PXP Hamiltonian and initial states being the z-direction computational bases, a x (t) ≡ 0. We plot the rotation dynamics of (a y (t), a z (t)) in Fig. 5(a), (b) with initial states being |Z 2 and σ x L/2 |Z 2 respectively. Note that the norm of a(t) is not a constant number due to the entanglement between different spins. We can readily observe that the rotation speeds of spins are not constant, which are determined by the orientations of neighbor spins.
In plot (a), for the initial state |Z 2 , all the spins rotate anticlockwise and the spin configuration has a two-site translational symmetry (indicated by red and blue colors). On the contrary, in plot (b) for the initial state σ x L/2 |Z 2 , we find some strange inverse rotation regions (spins rotate clockwise for a short period). Interestingly, the positions of inverse rotation regions match up with the contour of the linear light cone to a good accuracy. These strange inverse rotations are a direct result of the kinetic constraints: For the three middle spins of σ x L/2 |Z 2 , | · · · ↑↓↓↓↑ · · · , they are forbidden to simultaneously rotate to the | ↑ state, resulting in the retarded rotations of the j = 10 and j = 12 spins (see plot (c) for a zoom-in image). Moreover, this retarded rotation behavior propagates outwards like a row of dominoes, forming a series of inverse rotation regions indicated by the dashed circles. The physical picture is as follows: On the anti-ferromagnetic background |Z 2 , we locally create a domain wall in the middle of the spin chain | · · · ↑↓↓↓↑ · · · , the quasi-particle travels ballistically in the following quench dynamics [21], which is consistent with the local thermalization picture discussed in the previous section.
In plot (b), after the inverse rotation regions, the anticlockwise kinetically constrained rotations are no more affected, and the two-site translational symmetry is restored inside the light cone. Because of the retarded rotations, the spin orientations inside the light cone become periodically distinguishable for the |Z 2 and σ x L/2 |Z 2 cases, leading to the periodic oscillation pattern of Holevo information (compare Fig.  5(c) and Fig. 2(e)). The analyses above once again confirm the conclusion that: The unusual information revival dynamics inside the light cone are due to the local rather than global thermalization and robustness of scarred eigenstates against local perturbations. After the wavefront of thermalization effect passes, the kinetically constrained rotations are preserved, causing the persistent oscillations of Holevo information in-  The half-chain entanglement entropy S for each energy eigenstate. The scarred eigenstates form a special band in the top of plot (a). The scars together with two thermal eigenstates in the bulk are marked by labels consistent with those in Fig. 7. side the light cone.

Appendix C: Numerical Methods and More Results
In this section, we illustrate the numerical methods used in the main text and provide more numerical results of the OTOC and Holevo information dynamics in quantum manybody scarred systems. We numerically simulate the ZZand XZ-OTOC dynamics of the PXP model using the timesplitting matrix product operator (MPO) method [30,84] up to system size L = 41, with the maximum bond dimension 300 and a Trotter step dt = 0.05. Fig. 6 displays a pictorial illustration of the algorithm. Specifically, we rewrite the OTOC The half-chain entanglement entropy S for each energy eigenstate. The scarred eigenstates form a special band in the top of plot (a). The scars together with two thermal eigenstates in the bulk are marked by labels consistent with those in Fig. 9.
where |ψ(t/2) = e −iHt/2 |ψ . Usually in order to maintain the simulation accuracy, we need to increase the bond dimension when extending the evolution time t, because of the increase of entanglement entropy. Here by equally splitting the Hamiltonian evolution to each part of F ij (t), we are able to reduce the support of the scrambled operators V j (t). For a fixed evolution time, this algorithm leads to notable reduction of the required bond dimension compared with previous MPO algorithms [30,86]. We use the time-evolving block decimation (TEBD) algorithm [85,86] based on the matrix product state (MPS) ansatz to simulate the Holevo information dynamics for the PXP model up to system size L = 41, with the maximum bond dimension 100 and a Trotter step dt = 0.05. The relatively low entanglement entropy of scarred eigenstates [ Fig. 8(b) and 10(b)] makes the simulations for the |Z 2 initial state much more efficient, yet this advantage no more exists for the |0 case. All the MPO and MPS based numerical simulations are carried out with the ITensor library [113].
Numerical calculations for small system size (L ∼ 20) are performed with the exact diagonalization (ED) method in the constrained Hilbert space of the PXP model, for instance, Fig.  2(f) and Fig. 3.
In the main text, we have shown that the OTOC dynamics exhibit exotic periodic revivals inside the light cone for the initial state |Z 2 . In Fig. 7 and Fig. 9, we present the ZZ-OTOC dynamics of the PXP model for more general initial states in the open boundary condition (with boundary terms σ x 1 P 2 and P L−1 σ x L ) and periodic boundary condition, respectively. First the ZZ-OTOC dynamics for superposition states of |Z 2 and |Z 2 = i σ x i |Z 2 are displayed in Fig. 7(a), (b). Second, in Fig. 7(c)-(i), we take the energy eigenstates of the PXP Hamiltonian as initial states for the ZZ-OTOC calculations. The eigenstates are marked with the corresponding labels in Fig. 8. We observe that despite lower oscillation contrast compared to the |Z 2 case, the pattern of persistent and synchronized oscillations for ZZ-OTOCs still exists for the superposition states of |Z 2 and |Z 2 (a)-(b), scarred energy eigenstates (c)-(g), and also general superposition states of scarred eigenstates. However, the oscillation pattern is absent for typical thermal eigenstates (h)-(i). The numerical results indicate that the periodic revivals of OTOCs inside the light cone are not some fine-tuned results for the initial state |Z 2 , but a general phenomenon for a certain class of states within the non-thermal scarred subspace. The information scrambling dynamics for quantum many-body scarred systems are intrinsically different from the thermal or many-body localized systems. One additional remark is that we have calculated the ZZ-OTOC dynamics for all the scarred eigenstates marked in Fig. 8 and found that the periodic oscillation pattern always exists while the oscillation contrast fades away for a few scarred eigenstates near the ground state, which is probably due to the finite size effect.
In order to distinguish the scarred eigenstates and typical thermal eigenstates, we show the overlap between the |Z 2 state and each energy eigenstate |E n of the PXP Hamiltonian with the open boundary condition in Fig. 8(a), and the half-chain entanglement entropy S of energy eigenstates in Fig. 8(b) [12,13]. All the scarred eigenstates and two thermal eigenstates in the bulk are marked by labels consistent with those in Fig. 7. The half-chain entanglement entropy of two scarred eigenstates near E = 0 is relatively larger compared with other scarred eigenstates, which is probably due to the open boundary condition. In order to rule out the possible effects of boundary conditions, we display the corresponding results of the PXP model with the periodic boundary condition in Fig. 9 and Fig. 10. Compared with Fig. 7 and Fig. 8, these results do not show distinct differences despite different boundary conditions. We present more numerical simulations of ZZ-OTOC and Holevo information dynamics of the experimental Rydbergatom Hamiltonian (Eq. (5) and related H ± in the main text) in Fig. 11 to demonstrate their measurable signatures. In plots (a), (d), the experimental parameters are the same as those used in [10], where Ω = 2, U 1 = 24, U 2 = 0.38 and ∆ = 0. We observe that the periodic oscillations display a much lower contrast than that of the PXP model, which is induced by the invalidity of the condition U 1 Ω U 2 . Constrained by the geometry of 1D equally-spaced atoms, U 1 /U 2 ∼ 64 is fixed, such that Ω around √ U 1 U 2 can best fulfill the condition above. With this motivation, we show the results in plots (b), (e), with parameters Ω = 1.5, U 1 = 12, U 2 = 0.19 and ∆ = 0. 19. We observe that the oscillation pattern could be identified more clearly. As mentioned in the main text, we have further added a non-zero detuning ∆, in order to offset the influence induced by U 2 . Through the simple grid search optimization of the parameters Ω and ∆, in plots (c), (f), among several instances with relatively high oscillation contrast, we display the one with parameters Ω = 2, U 1 = 12, U 2 = 0.19 and ∆ = 0.38, which is the 2D view of Fig. 3 in the main text. The linear light cone contour and periodic oscillations inside the light cone can be readily observed.
In the main text, we have numerically calculated the XZ-OTOC dynamics of the PXP model. However, XZ-OTOCs can not be directly reduced to the observable average form like Eq. (6) in the main text, thus do not have a similar measurement scheme like the ZZ-OTOCs. Several other methods such as the ones based on randomized measurements [114] or classical shadow estimations [115] might be modified and adopted instead.