Fig. 2: Robust topological edge modes at up to infinite temperature.

a, Quantum circuit for implementing U(δt), which emulates a single-step evolution (288 ns) under the Hamiltonian in equation (1). The system is initialized in either the manifold {|Ψ₀⟩} (excitation number n = 0, corresponding to zero temperature), the manifold {|Ψ_e⟩} (n ≠ 0, finite temperature), or the product states \(| \bullet 00\ldots 0\bullet \rangle \) (effectively infinite temperature), and evolved with U(δt) for t cycles. Here, J_o, J_e and h_x are parameterized into the rotation angle θ around the x axis of the Bloch sphere [X(θ)]. V_xx is encoded in a combination of controlled-phase gates [CPhase(−2V_xx)] and Z phase gates [Z(V_xx)]. b, Measured time dynamics for the left edge operators in the homogeneous case (J_o = J_e = π/5). Black lines show the results of echo circuits, which estimate the decay caused by circuit errors (see Supplementary Information section 2G for more data and discussions). c, Measured site-resolved dynamics of normalized expectation value \(\overline{\langle {K}_{i}\rangle }\) for bulk stabilizers \(\{{\sigma }_{i-1}^{z}{\sigma }_{i}^{x}{\sigma }_{i+1}^{z}\}\) and edge operator \({\mathop{X}\limits^{ \sim }}_{{\rm{L}}}\) in the homogeneous case (J_o = J_e = π/5) near the left edge. The nearest excitations to the left edge are initialized at {Q₃, Q₅} (top) and {Q₅, Q₇} (bottom). d, Measured time dynamics of the left edge operators with fixed J_e = π/5 and varying J_o. Resonant processes lead to enhanced decay rates at J_o/J_e = 1 for \({\widetilde{Z}}_{{\rm{L}}}\) and J_o/J_e = 1, 2 for \({\widetilde{X}}_{{\rm{L}}}\). Error bars in b and d represent the standard deviation over five rounds of measurements, with each taking 10,000 shots. The time dynamics of the right-edge operators are shown in Extended Data Fig. 3. e, Spatial profile of the prethermal strong zero mode \({\varPsi }_{{\rm{L}}}^{z}\). The solid boxes denote theoretical predictions, with black frames highlighting the positive values and red frames highlighting the negative values. The coefficients are obtained by averaging the late-time dynamics over cycles from t = 25–40, with the sum of their squares normalized to unity.