Disorder-tunable entanglement at infinite temperature

Emerging quantum technologies hold the promise of unravelling difficult problems ranging from condensed matter to high-energy physics while, at the same time, motivating the search for unprecedented phenomena in their setting. Here, we use a custom-built superconducting qubit ladder to realize non-thermalizing states with rich entanglement structures in the middle of the energy spectrum. Despite effectively forming an “infinite” temperature ensemble, these states robustly encode quantum information far from equilibrium, as we demonstrate by measuring the fidelity and entanglement entropy in the quench dynamics of the ladder. Our approach harnesses the recently proposed type of non-ergodic behavior known as “rainbow scar,” which allows us to obtain analytically exact eigenfunctions whose ergodicity-breaking properties can be conveniently controlled by randomizing the couplings of the model without affecting their energy. The on-demand tunability of quantum correlations via disorder allows for in situ control over ergodicity breaking, and it provides a knob for designing exotic many-body states that defy thermalization.

The abundance of entanglement and other types of correlations in many-body systems make them an attractive resource for quantum information processing.Quantum coherence, however, is typically fragile: it rapidly deteriorates once the system is excited to a finite energy density above its ground state.This is a consequence of thermalization -the ultimate fate of generic systems comprising many interacting degrees of freedom [1][2][3].If thermalization breaks down, new types of dynamical behavior and phases of matter can emerge.For example, finely tuned one-dimensional systems [4] can evade thermalization due to their rich symmetry structure known as quantum integrability [5].On the other hand, in real materials disorder is ubiquitous and, if strong enough, it can strongly suppress thermalization by turning the system into an Anderson insulator [6] or its interacting cousin, the many-body localized (MBL) phase [7,8].
The ability to suppress thermalization while retaining a high degree of control over entanglement is key to technological applications.Integrable systems do not meet these requirements as they are restricted to one spatial dimension and require fine tuning of the system's parameters.In MBL systems, bulk excitations are localized, regardless of their energy density, which indeed may effectively protect the information stored in the degrees of freedom at the system's boundary [9,10].Nevertheless, entanglement in MBL states is typically bounded by "area-law" scaling [11].This limits their applications in quantum-enhanced metrology, which often rely on large multipartite entanglement [12].The latter entan-glement structures have recently been identified [13,14] in a class of systems known as quantum many-body scars (QMBS) [15][16][17].When QMBS systems are prepared in special initial states, their dynamics become trapped in a subspace that does not mix with the thermalizing bulk of the spectrum, leading to coherent time evolution of local observables [18][19][20][21][22][23].The observation of QMBS has triggered a flurry of theoretical efforts to understand and classify the general mechanisms of weak ergodicity breaking in isolated quantum systems [24][25][26][27][28][29][30][31][32][33].
In this work, we leverage the essence of QMBS phenomenon to create a new physical platform that supports a variety of entanglement structures that persist far from equilibrium and can be deterministically tuned by disorder.Our model is realized on a state-of-the-art superconducting qubit chip, in which the qubit-qubit coupling can be broadly tuned to encompass opposite coupling signs.The approach is inspired by the rainbow scar construction [34,35], which creates Bell pairs between qubits belonging to two halves of the system.We show that our model hosts several distinct families of QMBS states and entanglement structures.While the first family is a direct realization of the rainbow construction, the second family emerges from a hitherto unexplored mechanism: it is obtained by acting on the first family with the Hamiltonian of a single subsystem.By making the couplings spatially inhomogeneous, we can then turn these states into disordered QMBS states whose exact wave functions and entanglement structure can still be written down in the analytic form.Unlike their energies, the structure of these exact eigenstates can be explicitly modulated via the disorder profile, allowing to tune their properties.We experimentally observe the two types of entanglement via their distinct ergodicity-breaking signatures in quantum dynamics at late times.
The model and its symmetries.-Ourmodel is defined on a superconducting quantum processor with N = 2M qubits arranged in a ladder configuration, as shown in Fig. 1(a), with two horizontal rungs containing M qubits each.The coupling strength between a pair of nearest-neighbour transmon qubits can be capacitively tuned by a coupler, enabling broad ranges of [−8, 8] MHz and [−8, −2] MHz for the parallel and vertical couplings, respectively.Note that the unit 'MHz' has been divided by 2π throughout the entire text.We denote the states of each qubit by |•⟩ and |•⟩.The qubits belonging to the top row are described by ûα Pauli matrices (with α = x, y, z), while dα are Pauli matrices acting on the bottom-row qubits.The Hamiltonian can be written as where the top/bottom row and inter-row Hamiltonians, respectively, are given by Ĥσ=u,d = Here, the intralayer coupling amplitude J e,k and the frequency ω k can be site-dependent, allowing for the possibility of disorder, while inter-row coupling J a is assumed uniform.The rainbow construction mandates that the bottom row of qubits must have intra-row coupling amplitudes of the opposite sign.One can verify that the two Hamiltonians are related by the mirror transformation Ĥd = −M Ĥ⋆ u M † , where M simply maps , see Supplementary Information (SI) [36].Crucially, the mirror transformation forces the spectra of Ĥu and Ĥd to be identical but of opposite signs.
The model in Eq. ( 1) possesses a U(1) symmetry corresponding to the conservation of total magnetization along the z-direction and below, unless specified otherwise, we will restrict to its largest sector with zero magnetization.Interestingly, the exchange rules of the Hamiltonian give rise to an additional, more subtle, symmetry.Our Hamiltonian can exchange neighboring triplet states, T ≡ on each side in which the T and S are exchanged.Thus, Ĥ conserves the difference between the number of triplets and singlets, Q = M k=1 (n T,k − nS,k )(±1) l<k (n H,l +n D,l ) , multiplied by the phase factor that counts the number of doublons and holons to the left of a given site [36].The Q symmetry further splits the half-filling sector of the Hilbert space into M +1 disconnected sectors with quantum numbers −M, −M + 2, . . ., M − 2, M .Working in the largest Q sector (at zero magnetization), we checked the statistics of the energy level spacings of Ĥ using exact diagonalization, finding excellent agreement with the Wigner-Dyson ensemble and very small fluctuations between different disorder realizations [36], thus indicating that the model is quantum-chaotic.
Two families of rainbow scar entanglements.-Therainbow state |I⟩ = |TT . . .T⟩ is an eigenstate of our model in Eq. (1).Two different families of scars can be built from it using operators that commute with the halfsystem Hamiltonian Ĥu .To construct the first family, we use the operator Ẑ = M k=1 ûz k , which clearly commutes with Ĥu as the latter conserves z-magnetization.The powers of Ẑ are linearly independent up to ẐM , thus we can apply Ẑ up to M times.The Ẑ operator simply converts triplets into singlets and vice versa.The resulting states |M − n, n⟩ will be symmetric superpositions of all states with a fixed number n of triplets and M −n singlets.The scarred states of the first family are precisely these states, up to normalization: where n ranges between 0 and M −1 and |E M ⟩ ≡ |I⟩.The second equality illustrates that we can build |E n ⟩ recursively from |E n+1 ⟩ by making use of the projector P Q q on the sector of Q with an eigenvalue q.The projector P Q q is introduced for convenience as it simplifies the recursion, since acting with Ẑ on |E n ⟩ creates a superposition of |E n−1 ⟩ and |E n+1 ⟩.It can be verified that the state |E n ⟩ is an eigenstate of Ĥ with energy E n = J a (2n − M ) [36].While symmetry generators come to mind when we look for operators Ô that commute with Ĥu , we can build a different family of scarred states by using Ĥu itself as the generator.This gives us the second family of scars: with n = 1, 2, . . ., M − 1.The second term in the square bracket automatically orthogonalizes the states |E ′ n ⟩ with respect to the first family of scars.The projectors P Q q once again isolate E ′ n−1 from E ′ n+1 .Similar to the first type of scar, the second type of scar is also equidis-tant in energy, in fact occurring at the same energies E n = J a (2n − M ) [36].
It is instructive to contrast the two families of scars, Eq. ( 3) and (4).While both families occur at the same, regularly spaced energies throughout the spectrum, the total number of states in the second family is smaller by two than in the first family.Furthermore, there are stark differences in entanglement structures.The states belonging to the second family explicitly depend on disorder through their dependence on J e,k and ω k , unlike the first family.The states in the first family contain only singlets or triplets, with no doublons or holons, Fig. 1 The two families of scars are further identified by their entanglement entropy, S A =−trρ A log ρ A , where ρ A is the reduced density matrix of the subsystem A. The reduced density matrix ρ A =tr Āρ is obtained from the full density matrix ρ by tracing out the degrees of freedom belonging to the complement Ā of the subsystem A. The rainbow entanglement manifests as a striking difference in entropy depending on the type of bipartition that defines the subsystem A, and we will consider two types illustrated in Fig. 2(a).For the parallel cut between the two rows of the ladder, i.e., when A comprises qubits {u 1 , u 2 , . . ., u M }, the entanglement is large, as the bipartition cuts through an extensive number of Bell pairs.By contrast, for the bipartition perpendicular to the ladder, i.e., when , the entanglement is much lower.As seen in Fig. 2(b)-(c), this distinction is particularly striking for |E 0 ⟩ and |E M ⟩ states, which have nearly maximal entropy (i.e., scaling with the number of qubits) in the case of parallel bipartition but low entanglement (i.e., bounded by a constant) in the case of a perpendicular bipartition.
Entanglement dynamics.-Toobserve rainbow entanglements, we study global quench dynamics of the circuit, consisting of two contiguous rows with up to 8 qubits each.The structure of the first family of scarred eigenstates implies that they have high overlap with the product state |Π⟩ =        • .Figure 3(a) shows the dynamics of population imbalance, The population imbalance in the |Π⟩ state exhibits remarkable oscillations that persist up to time scales ∼1 µs.This is in contrast with typical thermalizing states, for which the population imbalance rapidly decays to zero by ∼50 ns.We note that it can be proven that our model in Eq. ( 1) exhibits perfect revivals of the imbalance from the |Π⟩ state [36].By contrast, a detailed characterization of the  The inset in (c) shows the entropy of a subsystem consisting of 2 qubits, sketched on the left.Purple dots with error bars stand for the average and standard deviation over 8 disorder realizations of J e,k , randomly selected from an interval ∈ [1, 3] MHz.Also, for reference, we show a typical thermalization case with an initial state of The lines show the results of numerical simulations for the same parameters, including the additional cross couplings Jx ≈ 0.3 MHz and nonlinearity of qubits η ≈ −175 MHz, present in the physical device [36].
experimental device shows that there are two extraneous couplings that are not present in the model [36].These are responsible for the slow population decay, observed in Fig. 3(a), as reproduced by the numerical simulations.The effect of these perturbations, however, is sufficiently weak such that clear signatures of the two families of scars can be observed and sharply distinguished.For example, a salient feature of the first family of scars is their insensitivity to the disorder in the Hamiltonian couplings.Thus, we expect the coherent dynamics from the |Π⟩ state to be unchanged when inhomogeneity is introduced in J e,k couplings.This signature is clearly confirmed by experimental observations in Fig. 3(a).
Furthermore, we utilized the quantum tomography technique and obtained the reduced density matrix of the subsystem consisting of qubits A = {k = 1, 2}, which gives us additional information about the dynam-ics beyond local observables.The subsystem fidelity, √ F A = tr ρ A (t)ρ A (0) ρ A (t), and entanglement entropy, S k=1,2 , are shown in Figs.3(b)-(c).For the initial state |Π⟩, the subsystem fidelity dynamics undergoes persistent revivals, implying that the initial information is restored many times, with a period of about 80 ns.Meanwhile, for generic initial product states, F k=1,2 quickly decays towards a value close to the inverse of the subsystem Hilbert space dimension, as shown in Fig. 3(b).The growth of entanglement entropy also shows a stark contrast between initial states.Compared to the thermal states, the |Π⟩ state exhibits a slow linear growth with superposed oscillations.The small oscillations are in correspondence with the peaks and valleys observed in the fidelity dynamics, with roughly half the period of the latter.Furthermore, we show the fidelity and entropy dynamics with different subsystems {u 1 , u 2 } and {u 1 , d 1 } in the inset of Fig. 3(c), confirming the rainbow entanglement structure previously sketched in Fig. 1(c).
To probe the second scar family, we require a more complicated initial state with an entangled doublonholon pair: • , which has predominant overlaps with the second family of scarred eigenstates [36].We emphasize that the state |ϕ L ⟩ is orthogonal to the first family of scars, as the latter do not contain any doublons or holons.To prepare the state |ϕ L ⟩, we use the circuit scheme in Fig. 4(a), which is composed of a few single-qubit and two-qubit gates.By utilizing high-precision tomography measurements, we then obtain the reduced density matrix of the subsystems {k = 1, 2} or {k = 2, 3}.The former one at t = 0 is visualized in Fig. 4(b), demonstrating that the entangled state |ϕ L ⟩ is successfully prepared.
To reveal the difference between the first and second family of scars, we focus on the subsystem A ′ ≡ {k = 2, 3}, whose fidelity F k=2,3 (t) and entanglement entropy S k=2,3 (t), are plotted in Figs.4(c)-(d).The subsystem fidelity partially reveals the scarred eigenstates and the Fourier transformation of F k=2,3 for the state |Π⟩ has an additional peak compared to the state |ϕ L ⟩.This difference is related to the fact that the first family of scarred eigenstates contains two more members compared to the second family in Fig. 2(b).Since only four qubits {k = 2, 3} are observed and the coherent information of the rest of the system is traced out, the dynamics of F k=2,3 (t) approximately reflects the subsystem itself.Thus, the Fourier spectrum of F k=2,3 (t) involves two peaks and one peak for state Π and ϕ L due to only two and one dimers, respectively [36].Furthermore, the choice of the subsystem is motivated by the fact that it leads to entropy ln 2 for the |ϕ L ⟩ initial state, while the entropy is still trivially zero for the |Π⟩ state.This distinction is verified in our experiment, as shown in Fig. 4(d).
Discussion and outlook.-Wehave realized a new platform that hosts multiple families of non-thermalizing states with distinct rainbow entanglement structures.Experimental distinction between the first and second family of scars.(a) Circuit diagram for generating the entangled initial state |ϕL⟩ used to probe the dynamics of the second scar family.Symbols "+","H", and "X" stand for CNOT, Hadamard, and Pauli-X gates, respectively.(b) Absolute values of the reduced density matrix elements ρ k=1,2 at t = 0 and the color of the elements means the phase.(c)-(d) Fidelity and entanglement entropy dynamics of a 4-qubit subsystem for initial states |Π⟩ and |ϕL⟩, which reflect the first and second family of scars, respectively.Insets of (c) show the Fourier spectra of the fidelity, which distinguish the second family (only one peak) from the first family (two peaks).The superconducting ladder contains N = 16 qubits with the same parameters as Fig. 3.
Our construction allows us to write down exact wave functions for these states, even in the presence of disorder.The existence and stability of scar states in disordered models have recently attracted much attention in theoretical studies [38][39][40][41][42].A unique aspect of our superconducting ladder circuit is that it allows us to tailor the entanglement of scarred states, while the explicit eigenstate dependence on the disorder profile controls the extent of ergodicity breaking.We observed signatures of rainbow entanglements by performing quantum state tomography of a many-body state of the ladder following the quench from special initial states, demonstrating robust revivals in the fidelity dynamics and slow growth of entanglement entropy far from equilibrium.
While in this work the disorder strength was assumed to be sufficiently weak such that the system overall remains chaotic, our versatile circuit allows direct access to strong ergodicity breaking regimes, where many-body lo-calization was recently proposed to give rise to new types of phenomena such as "inverted scarring" [43][44][45].More generally, our work bridges the gap between theoretical studies of QMBS, which place the emphasis on exact constructions of scarred eigenstates [16,17], and experimental realizations, e.g., in Rydberg atom arrays [18,19] or optical lattices [21], in which the scarred states are not known exactly (apart from a few exceptions [46]).In contrast, our model in Eq. ( 1) hosts exact scars, while its experimental implementation contains additional terms that act as a perturbation.While these perturbations were shown to be sufficiently weak in our device to allow unambiguous observation of scar signatures, it would be interesting to study in detail their effects on the stability of QMBS states in larger circuits that exceed the capability of classical simulations.
The flexibility of our construction stems from the fact that the scarred states are not generated by conventional symmetry action, but by the Hamiltonian describing one of the rungs of the ladder.For simplicity, we assumed the latter describes an integrable XY model (although the system overall is non-integrable).This was not essential, however, and our construction can be straightforwardly generalized to cases where the subsystem Hamiltonian is non-integrable [36].The key to implementing the construction was the broad tunability of the experimental device that allowed us to vary the coupling sign, in contrast with traditional multi-qubit superconducting systems [47].This tunability offers a new physical freedom that can be used for designing exotic many-body states that defy quantum thermalization and information scrambling, in particular the creation of multipartite-entangled states without the need for larger spins or complicated interactions [14].
inputs from other authors.All authors contributed to the experimental setup, discussions of the results and development of the manuscript.The experiment is performed on a flip-chip superconducting quantum processor hosting 2 × 20 frequencytunable transmon qubits in a ladder configuration.Each pair of adjacent qubits (both in rungs and side rails) are coupled by a tunable transmon coupler.All the control lines and readout resonators are located on a silicon substrate (bottom chip) and all the qubits and couplers are located on a sapphire substrate (top chip).The substrates are connected together with indium bumps as described in Ref. [48].Each qubit can be controlled by three control lines.The microwave (XY) line can control the state of the qubit, while the flux (DC/Z) lines can change the frequency of the qubit.The maximum resonance frequencies for the qubits are around 4.6 GHz and can be effectively tuned to below 4.0 GHz.The qubit can be measured using a capacitively coupled readout resonator, with the frequency of the resonator ranging from 6.5 GHz to 6.7 GHz.Each coupler is equipped with its own flux (DC/Z) lines to tune its frequency, with a designed maximum frequency of around 9 GHz.The total coupling strength of each pair of adjacent qubits is composed of direct coupling between them and indirect coupling through the coupler.The former is dominated by the direct capacitance between the qubits while the latter is determined by the capacitance between the qubit and the coupler.Both have been carefully designed so as to realize different tunable ranges in the rungs and side rails.Detailed information about the processor can be found in Tab.S1, including the idle frequencies, average single-qubit gate error, energy relaxation times, and dephasing times.Note that in the experiment, the qubits form a coupled system that is insensitive to each qubit's flux noise.Thus effective dephasing times are usually much longer than the Ramsey dephasing time T * 2 .For example, using the spin echo technique, the dephasing time is typically measured to be around 10 µs.

B. Experimental protocols
We use the quench protocol to observe the many-body scar on the processor.This includes three main steps: state preparation, interaction and measurement.In the experimental sequence, we first initialize all the qubits, Q i , in their ground state at their idle frequency ω j , the vertical couplers at around their sweet spots, and the horizontal couplers at around the frequency where the coupling strength between two nearest qubits are zero.
The initial states are prepared by alternating the single-qubit gates and the two-qubit controlled-π (CZ) gates.The preparation of the first family of scarred state |Π⟩, which is a product state, is realized by applying single-qubit XY rotations to every qubit.The preparation of the second family of scarred state |ϕ L ⟩ differs from the first family, as it is a "cat" state.We alternate singleand three two-qubit gates on the former four qubits, and apply XY rotations on the other qubits (as illustrated in Fig. 4 of the main text), with each coupler dynamically switched between nearly off and on to realize the CZ gate.CNOT gates in the experiment circuits are realized by CZ gate with Hadamard gate operated on the target qubit before and after.
In the step of interaction, we bias all qubits on resonance at the interaction frequency ω i .Meanwhile, the Table S1.Parameters of the processor: ω 0 j is the maximum frequency of Qj, known as the sweet spot; ω i j is the idle frequency of Qj where single-qubit XY rotations are applied and the average single-qubit gate error esq is measured via randomized benchmarking.All qubits are biased to ω I ≈ 4.375 GHz to activate the interaction, where the energy relaxation time T1,j and the Ramsey dephasing time T * 2,j of each qubit Qj are measured.Here, all frequencies have been divided by 2π.
Qubit ω 0 j (GHz) ω i j (GHz) esq (%) couplers are also tuned to turn on the net couplings between neighboring qubits.After waiting for an interaction time t, we finally tune all qubit frequencies away from the interacting regime in order to measure the relevant quantities at their readout frequencies.Repeating this process with varying t allows us to obtain the dynamics of the system.Also, using the technique of quantum state tomography, we can get the reduced density matrix of the subsystem A.
In the experiment, tuning the coupling strengths of all couplers and tuning all qubits on resonance are two very important steps, for which we have conducted careful calibration following the procedure below: Figure S1.Typical tunable effective couplings.The upper and lower panels stand for a pair of horizontal qubits (u4 − u5) and a pair of longitudinal qubits (u4 − d4), respectively.The left panels show the two adjacent qubits swapping dynamics as tuned by the coupler.The right panels are the absolute effective coupling strengths as a function of the coupler Z pulse amplitude, obtained by fitting the oscillations on the left panels.The horizontal coupling strength can be adjusted from positive to negative values, while the longitudinal coupling strength can be tuned in a negative range.
1. Coarse-tune coupling strengths: the coupling strength between each pair of neighboring qubits is coarsely calibrated as a function of the amplitude of the coupler Z pulse, as shown in Fig. S1.In this process, we excite one of the qubits and then place them at ω i for 200 ns, during which other qubits are placed 50 MHz above ω i , and other couplers are placed at around their maximum frequency.The coupling strength can be estimated by fitting the swapping dynamics.Following this process, we obtain the functions of all couplers.
2. Fine-tune coupling strengths: Keeping the frequencies of other qubits as above, we apply the corresponding Z pulse on other couplers according to the above results to roughly achieve the designed coupling strength.Then we slightly change the coupler Z pulse to fine-tune the coupling strength, validated by the swapping dynamics.This process is conducted for each pair of neighboring qubits.
3. Fine-tune frequencies of qubits: The frequencies of qubits may be affected by some factors such as Z pulse distortion, so we adopt the following strategies to minimize this uncertainty, during which all couplers are placed at their desired frequencies.We first apply π/2 pulse on one qubit, then place other qubits at a frequency above ω i and tune this qubit to ω i for a fixed delay, and finally measure its accumulated phase ϕ + = ∆ω + × t by tomographic operations.We also measure its accumulated phase ϕ − = ∆ω − × t when other qubits are placed 50 MHz below ω i .We then slightly change the qubit frequency to make ϕ + + ϕ − ∼ 0. This process is implemented for every qubit.

II. EFFECT OF CROSS COUPLING AND NONLINEARITY IN THE DEVICE
The Hamiltonian describing our experimental device can be written as Here, Ĥ denotes the Hamiltonian of Eq. ( 1) in the main text, which has been reformulated in terms of bosons, where the standard raising and lowering operators are given by û± = (û x ± iû y ) /2 (and, similarly, for d± ).
We consider a maximum of 2 photons per site.The last two terms represent the experimental imperfections due to the cross coupling J x between the diagonal qubits and nonlinearity η of the transmon qubit [22].The exact value of the J x couplings as measured in our device are listed in Table S2.These terms have been taken into consideration in the numerical simulations presented in Figs.3-4 of the main text.The existence of these two terms weakens the amplitude of revival dynamics of both families of scar states.However, the perturbed model Ĥexp still supports the two scar families and allows us to clearly distinguish between them.In Fig. S2 we show the experimental results of imbalance, subsystem fidelity, and entanglement entropy of the first four qubits from system size of N = 10 to 16, in good agreement with the numerical simulations.
While the results in Fig. S2 already suggest excellent convergence of local observable properties with system size, in Fig. S3 we extend the computation of the same observables in much larger systems N ≤ 40 using matrix-product state methods and time-dependent variational principle (TDVP), as implemented in TenPy libraries [49].The accuracy of this calculation is controlled by the bond dimension of the matrix product states, and the convergence was ensured by requiring that the relative error in F k=2,3 (t) and in S k=2,3 (t) observables are always below 10 −3 when comparing the two largest bond-dimensions used.Global quantities, such as manybody fidelity and bipartite entanglement entropy, were also monitored and showed good agreement between the different bond dimensions.The strength of the cross couplings J x beyond the 16 first qubits have been randomly drawn from a uniform distribution in [0.05, 0.45] and are then kept identical across all system sizes.Consistent with the results in Fig. 4 in the main text, the fidelity dynamics of subsystem A = {k = 2, 3} involves two frequencies for the first family of scar, one more than in the second family.The initial entanglement entropy is also ln 2 for the second family of scar, distinguishing it from the first scar family which has zero initial value for the entropy.The data in Fig. S3 reveals a remarkably fast convergence in system size.While the smallest system N = 8 deviates slightly from larger sizes, we can see that F 2,3 (t) and S 2,3 (t) are already fully converged across the experimentally accessible time interval of 500ns for system sizes N ≥ 16.

III. SUBSYSTEM FIDELITY FOR TWO FAMILIES OF SCARS
In general, the fidelity of a subsystem A for both pure and mixed initial states is given by The reduced density matrix is given by where The first scar family can be probed by studying the dynamics of |Π⟩, which is constrained to a hypercube subspace with a single particle in each dimer {u k , d k }.Then, the dynamics of the reduced density matrix simplifies to Thus, the reduced density matrix of the subsystem A for |Π⟩ is equivalent to an isolated system of size N A .For the subsystem A = {k = 1}, its fidelity is a cosine function with only one peak in its Fourier spectrum.As the subsystem size increases, the revival fidelity consists of For the initial state |ϕ L ⟩, the fidelity revivals are not exact even for the ideal Hamiltonian Ĥ.In this case, the reduced matrix cannot be simplified as ρ A (Π; t) above.However, to understand the fidelity dynamics it is helpful to consider a small system.For A = {k = 1} and {k = 1, 2}, the corresponding initial state are , which are entangled states.Their fidelity dynamics have no revivals.If the subsystem A involves the third dimers, i.e.A = {k = 1, 2, 3}, the initial state Its fidelity dynamics have onefrequency revival dynamics, as shown in Fig. S5.
Furthermore, we emphasize that the experimental perturbations do not affect our distinction between states Π and ϕ L .In Figs.S4 and S5, we numerically compute the fidelity dynamics for different subsystem sizes based on the experimental Hamiltonian.The signatures that distinguish the first and second families of scars can be clearly observed.

IV. HAMILTONIAN IN THE DIMER BASIS
The Hamiltonian from the main text can be re-written as The ûα and dα are Pauli matrices acting on the top and bottom rows, respectively.The state • denotes an upspin and • denotes a down-spin, such that û The effect of applying the Hamiltonian on any state in the computational (Z) basis is straightforward.However, its effect on the dimer basis is not obvious and we explicitly work it out in this section.To achieve this, it will be convenient to use the above decomposition of the Hamiltonian as a sum of terms that either act on a single dimer or on two dimers at once.The former corresponds to diagonal terms and hopping perpendicular to the ladder (denoted by ĥ⊥ k ), while the latter encompasses hopping parallel to the ladder and is denoted by ĥ∥ k,k+1 .In the main text, we introduced the convenient local dimer basis, spanned by doublon |D⟩, holon |H⟩, triplet |T⟩ and singlet |S⟩ states.Those were given by: It is straightforward to see that every state in the dimer basis is an eigenstate of ĥ⊥ k with The action of ĥ|| k,k+1 is more complicated and requires looking at two neighboring dimers.It is straightforward (S9) The missing combinations can be obtained by flipping the two sites under consideration on the right-and left-hand sides.Importantly, this leads to the following combinations of states being annihilated by ĥ∥ k,k+1 : (S10)

V. SYMMETRIES OF THE MODEL AND STATISTICS OF ENERGY LEVELS
Our model, defined in Eq. (S5), conserves the quantity Q that was introduced in the main text: where Tk , Ŝk , Dk and Ĥk are projectors on the respective dimer state at the given site k.Here we outline the proof of this statement and provide physical intuition behind the conservation of this quantity.We start by considering the action of the Hamiltonian in the dimer basis, Eq. (S7).The relevant dynamical terms are summarized in Fig. S6(a).Consider a simple configuration like TSSSTS.The Hamiltonian rules dictate that we can only exchange TS or ST into HD or DH.Therefore, a natural guess for the conserved quantity would be the difference between the number of triplets and the number of singlets: T − Ŝ.This works perfectly until we start having configurations with neighboring dimers like SH.In that case, we can turn it into HT and thus change T − Ŝ.However, for that to happen it means that we first created HD and that the dimer that was switched between S and T must be surrounded by H and D. Effectively, the Hamiltonian creates domains inside which all T and S are "exchanged".Thus, we can keep track of how many times such exchanges have occurred for a given dimer by counting the number of doublons or holons to its left.Fig. S6(b) shows an example of this process.
While we focused on the half-filling sector in the main text, we note that Q is a symmetry at any filling.The interplay of Q and filling leads to a large number of disconnected sectors.Some of them are of small dimension and similar to the ones studied in Ref. [50].Nonetheless, for large sectors we find good agreement with the eigenstate thermalization hypothesis (ETH) predictions, as shown in Fig. S7.Finally, beyond the symmetry Q, at half-filling the Hamiltonian anticommutes with where Pd↔u is operator exchanging the top and bottom rows.Ĉ is related to a particle-hole transformation along with a swap between the top and bottom row and an additional phase.In the dimer picture, Ĉ simply switches T and S as well as D and H, with a (−1) phase for every D present.Note that, as we are at half filling, there are always the same number of D and H, so Ĉ = Ĉ † and Ĉ2 = 1.As it also anti-commutes with Q, Ĉ exchanges the sectors with Q = q and Q = −q.For M -even, we have a sector with q = 0; in that sector the spectrum is symmetric around E = 0 because of Ĉ, however, this is not the case in other sectors.

VI. STRUCTURE OF SCARRED EIGENSTATES
The rainbow scar construction [34,35] is based on two subsystems (labeled "1" and "2"), with respective Hamiltonians obeying the relation As the spectra of two subsystems are identical up to a minus sign, the composite system has a large zero-energy subspace, spanned by pairs of eigenstates with energies E and −E.Provided that M maps each site in subsystem 1 to a single site in subsystem 2, as in Eq. (S13) (i.e., with no global term mixing all sites), then the zero-energy subspace contains a rainbow state, where D 1 denotes the number of states in the subsystem 1 and σ j denotes the state of site j in it.The state |I⟩ is simply a tensor product of Bell pairs, and it has maximal entanglement between the two subsystems.Importantly, |I⟩ is an eigenstate of the full system, independently of the microscopic details of Ĥ1 [34].One can then make the two subsystems interact by adding a term Ĥint .If it is chosen such that it has |I⟩ as an eigenstate, |I⟩ will the be a scarred eigenstate of the full system while other states generically become ergodic.
In our case, subsystem 1 is the top row of the ladder while subsystem 2 is the bottom row.The transformation relating their Hamiltonians is where Pd↔u is the operator exchanging the top and bottom row and the term between parentheses is a particlehole exchange.The rainbow state can then be written succinctly as which reveals the structure of Bell pairs formed between the two rungs of the ladder.Upon addition of Ĥint , the rainbow state |I⟩ remains an exact eigenstate with energy M J a , thus becoming a rainbow scar of the full model.The previous construction is not limited to a single state.We can generate additional rainbow scar states by acting on |I⟩ with an operator Ô ⊗1 2 , where Ô commutes with the Hamiltonian Ĥ1 [34].The resulting state, by construction, still belongs to the zero-energy subspace of the combination of Ĥ1 and Ĥ2 .Provided that the resulting state is also an eigenstate of Ĥint , it is then a scarred eigenstate of the full system.Two different families of scars in our model in Eq. ( 1) can be built this way by choosing different operators Ô.Moreover, it allows to get simple product states with overlap only on these special eigenstates.

A. First scar family
To construct the first scar family, we use the operator Ẑ = M k=1 ûz k .As explained in the main text, this operator commutes with Ĥ1 and it can be applied up to M times, changing triplets into singlets and vice versa.The resulting states, |E n ⟩, are members of the first type of scar family.They represent a symmetric superposition of all configurations with a fixed number of singlets and triplets, and it is straightforward to verify that they are eigenstates of the model with energy The underlying structure of the first family of scarred states |E n ⟩ is the restricted spectrum generating algebra [29].This is seen by noting that they can be built using the raising operator Ŝ+ = M k=1 |T k ⟩ ⟨S k |.Thus, we can equivalently express the first type of scar states as where N is a normalization constant.We recognize that these scarred eigenstates are built in a similar way to previously known examples, e.g., in the spin-1 XY model [28].In Sec.X we compute their entanglement entropy for the bipartition perpendicular to the ladder, and we find that the eigenstate at zero energy has S 1,⊥ = 0.5 + 0.5 log(πM/8), in the limit M → ∞.We note that the su(2) algebraic structure of these states implies they have extensive multipartite entanglement [13,14].

B. Second scar family
The second scar family can be built recursively from the first family, where the role of the operator Ô in the rainbow construction is played by the subsystem Hamiltonian Ĥ1 .Specifically, the second type of scar states are given by raising from |E n−1 ⟩ according to or, equivalently, by lowering from |E n+1 ⟩ as In these expressions, P Q q is a projector on the sector of Q with eigenvalue q.
In order to write the states |E ′ n ⟩ in a more explicit form, let us define ensembles of sites Λ = {1, 2, . . .M }, and |a, b⟩ A the symmetric superposition of all states with a singlets and b triplets on sites in the ensemble A. More details about the construction of these states and the proof that they are eigenstates of the model is given in Secs.VII-VIII, where we also derive the normalization factors N n .There, we prove that |E ′ n ⟩ have the same energies as the first family of scar states, i.e., E ′ n = J a (2n − M ).
To the best of our knowledge, there is no su(2) algebraic construction for the second family of scarred states |E ′ n ⟩, despite their equal energy spacing.Indeed, the latter are generated by acting with an operator on |E n ⟩, and not on E ′ n−1 or E ′ n+1 .Moreover, the second type of scarred eigenstates depend sensitively on the parameters J e,k and ω k of the model.Therefore, the entanglement of these states is not fixed, like in the first family of scars, as we discuss in detail in Sec.X.For a cut perpendicular to the ladder, the second type of scarred states generally possess entanglement entropy that scales logarithmically with subsystem size.While we are unable to rule out the possibility of volume-law entanglement scaling, our analytical and numerical results strongly suggest that, in the large M limit, the second type of scarred state with zero energy obeys S ⊥,1 < S ⊥,2 < S ⊥,1 + log 4 [36].This suggests it should always be possible to find a low-entangled state which has overlap only on the scarred stats of the second family, regardless of the values of parameters.By tuning the hopping J e,k one can then change which state has this anomalous overlap, as we will demonstrate in Sec.IX.

VII. BUILDING SCARRED STATES OF THE SECOND FAMILY
In this section we derive the exact form of the scarred states of the second family.The formal proof that these states are eigenstates of the model is delegated to the subsequent section.
In order to reveal the structure of the second family of scarred states, let us first derive the action of Ĥ1 − ( k ω k /M ) Ẑ in the dimer basis.First we define ĥ∥,1 From there it is straightforward to see that where ω k = ω k − ( i ω i /M ).As we will only apply this operator to scarred states of the first family, which contain no |D⟩ or |H⟩, we can ignore any configurations containing them.The action of ĥ∥,1 k,k+1 and u z k on dimers is then From there, we immediately see that To represent the symmetric superposition, we recall the notation introduced to describe scarred states of the first family: This state is a fully symmetric superposition of all configurations on L sites with L − n singlets and n triplets.
For brevity, we will first work out what happens when we apply ĥ∥,1 1,2 (the same is true for any ĥ∥,1 k,k+1 ): (S26) where we used the condition (S24) to cancel the contribution of the first term.Applying the projector P Q M −2n singles out one of the terms: (S27) Next, we look at the action of ω 1 ûz 1 : (S28) Applying the projector in this case gives The result is similar if we act on another site k, with a prefactor ω k and a triplet on that site.Ultimately, we end up with a collection of all states with n singlets and n − 1 triplets, but each of them has a prefactor that depends on the location of the triplets.Let us introduce the operator Tk that gives 1 if this site is a triplet and 0 otherwise.We can then write To write down the scarred states of the second family, we now simply need to gather the terms from Eqs. (S27) and (S30).We also remove the overall minus sign and normalize the state.Let us introduce the ensembles of sites They represent, respectively, all sites, sites k and k+1, and all sites except Applying the Hamiltonian to this state and after some algebra, we obtain (S42) For |ψ⟩ to be an eigenstate, the coefficients must obey In general, we have 5 unknowns but only 3 equations, leaving room for two non-trivial solutions.The first option is to set all α j to be equal and all β j to zero: This corresponds to the scarred state of the first family, |E 1 ⟩.The other solution is given by the scarred states of the second family that obey (up to a normalization factor) Furthermore, we can verify that the two families of scars are orthogonal as their overlap is given by This completes the proof for the special case of M = 3.However, generalizing this to an arbitrary M is now straightforward.The general state is Now we have M different α j and M − 1 different β j , so 2M −1 unknowns in total.For each block of 2×2 sites we get an equation similar to Eqs. (S43) and (S44).Thus, for j = 1 to M − 1 we have For each rectangular block of 2 × 3 sites we get an equation similar to Eq. (S45).Hence, for j = 1 to M − 2 we have J e,j+1 β j − J e,j β j + 1 = 0. (S51) For a chain with 2M sites, this yields (M −1)+(M −2) = 2M − 3 equations.So for any system size we always get at least two scarred states with energy −(M − 2)J a .It is easy to check that the two solutions given in Eqs.(S46) and (S47) are still valid.

C. Other values n
While we treated the case n = 1 on its own to provide a simple example, the recipe is exactly the same for general n (except for n = 0 and n = M that were already proven).We use the same Ansatz in which we only consider states with one DH + HD pair in a background of n − 1 triplets and M − n − 1 singlets and states with n triplets and M − n singlets.
We first restrict our investigation to an arbitrary location k, k + 1 for the hole-doublon pair.
where X and Y denote either S or T. For ĥ⊥ , the contribution on all other sites except k and k + 1 is diagonal and equal to J a (2n − M ).Therefore, to prove that states are eigenstates with this energy, we need to prove that the action of the rest of the Hamiltonian annihilates the state.For ĥ∥ k,k+1 we only have to care about the action on sites k − 1 to k + 2. Indeed, the rest of the state is composed of triplets and singlets.For any other pair, if it is SS or TT, then it is annihilated by the action of the Hamiltonian.If, instead, it is TS, then there exists another state in the superposition, with the same weight, that has ST instead.Therefore, their superposition is also annihilated by the Hamiltonian.(S53) and so We get a unique equation for every of the M − 1 pair of sites and for every of the M −2 n−1 possible background configurations, where n = 1 to M − 1 is the index of the scarred states.Now we still need to look at the effect of ĥ∥ on XD, XH (as well as HY and DY).For that, we need to also consider the DH and HD pair placed one site to the left: (S55) Applying the Hamiltonian to these states leads to where X = T if X = S and X = S is X = T. Hence we get an equation for j = 1 to M −2.Thus, for the β coefficients, we always have M − 1 unknown and M − 2 equations.Furthermore, these equations take the form of Eq. (S57) and are identical for any value of n.Once again we recognize that setting all α equal and all β to 0 is a valid solution, and so the first family of scarred state is indeed an eigenstate.As for the n = 1, we also recognize that Eq. (S57) admits β j = J e,j /2 as a solution.
For the α, if we add a contribution of ω j for each site that has a triplet on site j as in the second family of scars, we recognise that α ST and α T S have the same contributions outside of site k and k + 1.It is then easy to see that α ST − α T S = ω j+1 − ω j and that Eq. (S54) is satisfied, concluding our proof.

IX. DYNAMICAL SIGNATURES OF TWO SCAR FAMILIES
Both families of scarred states are evenly spaced in energy with spacing 2J a , hence they can detected by persistent revivals following the quench from a suitably chosen initial state.In this section, we identify such initial states for the two scar families, and prove the existence of revivals by computing the overlap of initial states with scarred eigenstates derived in previous sections.We illustrate how these results can be used to enhance quantum revivals by modulating the coupling between the qubits.

A. Reviving initial state for the first scar family
To probe the first scar family, in the main text we used the state |Π⟩ given by: This state must undergo perfect revivals due to the exact su(2) algebraic structure, as we now explain.From the raising operator Ŝ+ = k |T k ⟩ ⟨S k |, we can infer Ŝx ∝ Ŝ+ + ( Ŝ+ ) † or, equivalently, The corresponding Ŝz operator is given by Ŝz = , which is easily seen to have the same action as Ĥ in this subspace.Hence, if we prepare the system in the lowest-weight state of Ŝx , it will undergo perfect precession around an effective field in the z-direction.It is easy to see that our state |Π⟩ in indeed the lowest weight state of Eq. (S59).The precession that starts out in |Π⟩ state will result in perfect state transfer to the highest-weight state • .The existence of revivals from the |Π⟩ initial state can be shown more explicitly by computing the overlap |⟨Π|E n ⟩| 2 with eigenstates of the first scar family.First, it will be important to observe that we obtain a simple product state when all possible combinations of singlets and triplets are summed up: (S60) where in the second equality we made use of The same procedure can applied to the |Π⟩ state by not- • .This means that there is now a factor of −1 for each singlet present, such that Thus, the |Π⟩ and |Π ′ ⟩ states only have overlap on the first family of scarred eigenstates.As the latter are regularly spaced in energy, the dynamics initialized in |Π⟩ or |Π ′ ⟩ exhibits perfect revivals.

B. Reviving initial state for the second scar family
For scarred states of the second family, the identification of the reviving initial state |ϕ L ⟩ is more subtle as the eigenstates now depend on disorder realization and there is no apparent algebraic structure.However, we can once again appeal to the fact that the sum of all singlet and triplet configurations can be written as a simple state in the Fock basis.For any n, |E ′ n ⟩ contain terms with N n the normalization factor given in Eq. (S38).Recall that Λ k denotes sites k and k + 1 and Λ k denotes all sites except k and k + 1.As a consequence, summing all (S64) In order to simplify this state and to make the orthogonality with |E n ⟩ obvious, we will consider the initial state with the normalization factor Z = 2 M −1 k=1 J 2 e,k .As the |ϕ ′ J ⟩ state only has overlap with states containing one doublon and one hole, it has zero overlap with the scarred states of the first family that have neither of those.As such, any non-trivial dynamics after a quench from this state must come from scarred states of the second family.We can compute this overlap exactly: which leads to It is important to notice here that it is not the ω k that enter this equation but their counterparts ω k , with the mean removed.So if we draw the J e,k and ω k from the same distribution [∆ − δ, ∆ + δ], if ∆ ≫ δ then we will have that J 2 e,k ≫ ω 2 k as the former is at the scale of ∆ 2 but the latter at the scale of δ 2 .In that case the state |ϕ ′ J ⟩ will have total overlap of order 1 + 2M δ 2 (M −1)∆ 2 −1 on the QMBSs of the second family, making these states the only relevant ones.
Similarly, we can also define (S68) This state obeys the same |⟨ϕ J |E ′ n ⟩| 2 as Eq.(S66), but with If the J e only have a small amount of disorder, we can also look at the homogeneous state Due to its homogeneous weights, it requires less fine tuning to prepare.It is well suited for the situation considered before, where [∆ − δ, ∆ + δ] with ∆ ≫ δ.
Finally, for experimental implementations, it is preferable to use an initial state that requires fewest gates to prepare.Because of this, in the main text we have considered the simpler cousin |ϕ L ⟩ of the state |ϕ⟩: For this state, the overlap with the scarred states of the second family is It is important to note that |ϕ ′ J ⟩, |ϕ J ⟩, |ϕ⟩ and |ϕ L ⟩ only have overlap with states with exactly one hole and one doublon.As a consequence they are completely orthogonal to scarred states of the first family which only have singlets and triplets.

C. Tunability of revivals
To confirm our analysis above, we have numerically computed the overlap of states |Π⟩, |ϕ J ⟩ and |ϕ L ⟩ with eigenstates of the model in Eq. (1) in Fig. S8.All the states have predominant support on scarred eigenstates, either of the first family (|Π⟩ state) or the second family (|ϕ J ⟩ and |ϕ L ⟩ states).We emphasize that, as scarred states of the first family have no doublons or holons, they are exactly orthogonal to |ϕ J ⟩ and |ϕ L ⟩. Thus, persistent revivals from the latter states provide unambiguous evidence for the second type of QMBS.Indeed, Fig. S9(a) shows that the initial states |Π⟩, |ϕ J ⟩ and |ϕ L ⟩ lead to revivals of the wave function and slow growth of entanglement entropy when compared to random Fock basis states -clear signatures of scarring.The rainbow nature of |E n ⟩ and |E ′ n ⟩ is also apparent in the dynamics, as the growth of entropy shows a stark difference between two different cuts, shown in Figs.S9(b)-(c).
While the revival fidelity from the |ϕ L ⟩ initial state in Fig. S9 is not particularly high, we can leverage the tunability of the second family of scarred states to enhance it.Indeed, from Eq. (S72) it directly follows that the projection of |ϕ L ⟩ state on the set of scarred eigenstates |E ′ n ⟩ is given by Thus, in the limit of J e,1 →∞, we recover perfect revivals.
. Improving the revivals by modulating the hopping strength on the first site according to Je,1 = J 0 e,1 + ∆1.Monotonic increase of revival peaks can be seen as ∆1 is increased.Data is for N = 16, Ja = 4, J e,k ∈ [4, 4.5], ω k ∈ [0.5, 1.5] drawn from a uniform distribution.
To illustrate the effect of J e,1 on the dynamics, in Fig. S10 we compute the fidelity dynamics as J e,1 is modulated by an amount ∆ 1 ∈ [0, 8].Setting ∆ 1 = 4 -which amounts to approximately doubling J e,1 -leads to a two-fold increase of the first revival peak, with the increase being even more pronounced at the second revival peak.This shows that the tunability of the initial state can be obtained with experimentally realistic values of the parameters for which the model remains chaotic.

X. ENTANGLEMENT ENTROPY OF SCARRED EIGENSTATES
Due to the su(2) algebra, the first family of scarred eigenstates have the structure of angular momentum eigenstates.As a result, they have entanglement entropy scaling as ∝ log M .In fact, it is straightforward to compute their entanglement entropy analytically.Let us assume that M is even and equal to 2R.For simplicity, we will concentrate on the state with n = R which has exactly zero energy, as it has the highest entanglement entropy among all scarred states.It is easy to decompose the state |E R ⟩ as where |ψ 1,k ⟩, |ψ 2,k ⟩ are the normalized versions of |R − k, k⟩ and |k, R − k⟩, respectively.From the last expression, we recognize the prefactors as the Schmidt coefficients.Therefore, the entanglement spectrum has R + 1 = M/2 + 1 non-zero values with for k = 0, 1, . . ., R. In the large-M limit, one can perform a saddle-point approximation to arrive at the result S 1,⊥ = 0.5 + 0.5 log(πM/8), demonstrating the logarithmic scaling with system size.
For the second family of scars, the computation is more arduous as the entanglement entropy depends on the disorder realization.Here we provide an exact computation for a few extremal cases.While we do not have proof that these are the cases with maximum and minimum entanglement entropy, they match with the results of our numerical optimizations.First, we can focus on the |J e,k | ≫ |ω k | case, in which we assume that the ω k are negligible.Let us also set all J e,k = 1 equal to 1, except for the middle one which we set to J e,R = J e .We can Consequently, the contribution to the entanglement spectrum in the first two sums are identical and given by Similarly, the third and fourth sums have the same coefficients given by The maximal entanglement entropy is obtained for J e = J ⋆ e and scales as √ M .In that case, we find numerically that S max 2,⊥ = S 1,⊥ + log 4 in the large-M limit.It is easy to understand how this additive factor can appear.For the first family of scarred states, we have R + 1 nonzero values in the entanglement spectrum, while in this case we have 4R − 2. For R very large, this is a fourfold increase of the number of non-zero values.As they have a similar distribution, this leads to a simple additive factor of 4 due to the log involved in the calculation.
The case with minimal entanglement entropy is in the limit of a single J e,k (not the middle) being much larger than all other ones.For simplicity, let us consider J e,k = δ 1,k .Then the state can be decomposed as (S83) In the limit of large M , we recover the same result as for scarred state of the first family S min 2,⊥ = S 1,⊥ = 0.5 + 0.5 log(πM/8).This can, once again, be understood simply from the number of nonzero values in the entanglement spectrum as their distribution is similar in both cases.As they have, respectively, R − 1 and R + 1 such values, they become identical at leading order in the large M limit.
For all three cases, we can compute the entanglement entropy efficiently from the analytical form of the entanglement spectrum for systems with hundreds of sites.Fig. S11 displays this along with the expected large-M behavior.We find very good agreement between them already at M ≈ 20, with the difference between them decreasing as 1/M .While we do not have a proof that the cases treated are the true maximum and minimum of S 2,⊥ we now illustrate numerically that they provide excellent bounds in This shows that high entanglement entropy can be obtained without getting close to a fine-tuned integrable point where all parameters are identical.
the large N limit.Fig. S12 shows the entanglement entropy obtained for multiple random realizations with different ranges of parameters, and for numerical minimization and maximization over all J e,k and ω k parameters.
For the minimum we find exact agreement between our analytical and numerical results.For the maximum case, we find that in smaller systems, realizations with high disorder can have higher entanglement entropy than our analytical Ansatz.Nonetheless, the difference between the numerical and analytical maxima converge as N gets larger.This is confirmed by looking at the entanglement spectrum.In all states obtained by numerical maximization of entropy, there are exactly N = 4R nonzero values in it.This precludes them from being volume-law states.
Asymptotically, this will also be equivalent to the 4R − 2 nonzero values in our analytical Ansatz, up to 1/N corrections.For these reasons, we believe that the scarred states of the second family cannot be volume-law entangled and that our analytical Ansatz provides an upperbound in the large N limit.The maxima of entanglement entropy obtained numerically are also useful to show that it cannot only be reached by fine-tuned cases where many parameters are identical.As these case might be integrable or have additional symmetries, they are usually not chaotic.Meanwhile, the numerical maxima do not have degenerate parameters and generically show characteristic of ergodic systems.This can be seen for example in the N = 12 case, for which the entanglement entropy of eigenstates is plotted in Fig. S12.The concentration of points around an arc is typical of systems obeying the ETH.

XI. GENERAL SYMMETRY SECTORS
We have studied different sectors of Q symmetry at non-zero magnetization, i.e., at general filling factors in the fermion representation of the model.Fig. S13 shows the entanglement entropy of eigenstates at half filling (zero magnetization).This is similar to the data presented in the main text, but we plot individual Q sectors for increased clarity.The scarred states are clearly visible, while the rest of the eigenstates show a narrow arc-like distribution, typical of a chaotic system.By contrast, Fig. S14 shows a sector away from half filling, where the same arc is visible.However, in this case no scarred state is present.Indeed, scarred eigenstates only exist at half-filling, and this can be directly understood from the rainbow scar construction.As the mirror transformation between the two subsystems M = ,.

XII. GENERALIZATION TO A NON-INTEGRABLE SUBSYSTEM HAMILTONIAN
In this section we demonstrate that our results can be generalized to the case in which the subsystem Hamiltonian Ĥ1 defines a non-integrable model.The model presented below is a disordered XY chain with nearestneighbor and next-nearest-neighbor couplings along the x and y directions.We use the same notation as in the main text and define Ĥ = (S84) This Hamiltonian also naturally decomposes into three parts: Ĥ1 (first two lines) acting on the top row of the ladder, Ĥ2 (third and fourth line) acting on the bottom

Figure 1 .
Figure 1.(a) Micrograph of the superconducting quantum processor in a ladder configuration.The tunable couplings, Ja and Je, between nearest-neighbor qubits, belonging to the same or opposite rows, are indicated.Blue and red curves illustrate the disordered coupling strengths ±J e,k , carrying opposite signs in the two rows.(b-c) Schematic representation of entanglement structure for a thermalizing state, the first family, and second family of scars, respectively.The shaded blue region in (b) indicates large entanglement between all qubits in the thermalizing case.In (c)-(d), dark curves depict the Bell pair entanglement of neighboring qubits, with dimers locally forming |T⟩ states.The tetramer configurations denote doublon-holon entanglement (|DH⟩ + |HD⟩) / √ 2, characteristic of the second scar family.

Figure 2 .
Figure 2. (a) Schematic of the ladder (M = 8) with dashed lines indicating two types of bipartitions.The parallel cut splits all entangled pairs in the rainbow state, while the perpendicular cut does not affect them.(b,c) Bipartite entropy of eigenstates in different Q-symmetry sectors for the two cuts shown in (a).The two families of scarred states are highlighted by red circles (first family) and blue diamonds (second family).Their rainbow nature is revealed by the fact that they move between minimal and maximal entropy, depending on the cut.The upper bounds of bipartite entropy in (b) and (c) are given by the Page entropy (2M ln 2 − 1)/2 [37] and the maximum subsystem entropy M ln 2, respectively.Colored cross sections represent different sectors labeled by the values on the Q axis.Data is obtained by exact diagonalization with N = 18 qubits, Ja = 4, and J e,k ∈ [4, 4.5], ω k ∈ [0.5, 1.5] drawn from a uniform distribution.
(c).By contrast, the second family has overlap with states involving a symmetric superposition on a single doublon-holon pair |. ..DH. ..⟩ + |. ..HD. ..⟩ with weight J e,k /2 on top of a background of M −n−1 singlets and n−1 triplets, Fig. 1(d).They also have overlap with all states with M −n singlets and n triplets, with prefactors depending on the ω k and on the location of the triplets.The dependence of |E ′ n ⟩ states on ω k and J e,k allows them to tune their properties, such as entanglement entropy or overlap with special initial states, as demonstrated below.

Figure 3 .
Figure 3. Experimental observation of the dynamical signatures of the first scar family.The ladder with N = 10 qubits is initialized in the state |Π⟩ and the couplings are set to values J e,k = Je = Ja/3 = −2 MHz.The plots show the measurements of (a) population imbalance, (b) the 4-qubit fidelity, and (c) the 4-qubit entanglement entropy, specified in the text.The inset in (c) shows the entropy of a subsystem consisting of 2 qubits, sketched on the left.Purple dots with error bars stand for the average and standard deviation over 8 disorder realizations of J e,k , randomly selected from an interval ∈ [1, 3] MHz.Also, for reference, we show a typical thermalization case with an initial state of •

F
Figure 4.Experimental distinction between the first and second family of scars.(a) Circuit diagram for generating the entangled initial state |ϕL⟩ used to probe the dynamics of the second scar family.Symbols "+","H", and "X" stand for CNOT, Hadamard, and Pauli-X gates, respectively.(b) Absolute values of the reduced density matrix elements ρ k=1,2 at t = 0 and the color of the elements means the phase.(c)-(d) Fidelity and entanglement entropy dynamics of a 4-qubit subsystem for initial states |Π⟩ and |ϕL⟩, which reflect the first and second family of scars, respectively.Insets of (c) show the Fourier spectra of the fidelity, which distinguish the second family (only one peak) from the first family (two peaks).The superconducting ladder contains N = 16 qubits with the same parameters as Fig.3.

D
Amplitude of fast Z pulse on coupler (a.u.)

Figure S2 .
Figure S2.Dynamics of imbalance (a,b), subsystem fidelity(c,d), and entanglement entropy (e,f) for simulations and experiments, respectively, with different system sizes ranging from N = 10 to 16.Data is for the first type of scar, with the same parameters as those used in Fig.3of the main text.

Figure S3 .
Figure S3.Influence of experimental perturbations.Fidelity and entanglement entropy dynamics of the subsystem {k = 2, 3} obtained by the numerical simulations of the ideal Hamiltonian Ĥ in the main text (panels a,b,e and f) and the Hamiltonian Ĥexp, which contains the perturbations in Eq. (S1) (panels c,d,g and h).Data for N ≥ 24 is obtained using TDVP.The coupling parameters are identical to those in Fig. 4 in the main text and are Ja = −6 MHz, J e,k = −2 MHz, η = 175 MHz and Jx ∈ [0.05, 0.45] MHz.

Figure S4 .Figure S5 .
Figure S4.Subsystem fidelity dynamics of the initial state |Π⟩ for Hamiltonians Ĥ in the main text and Ĥexp in Eq. (S1) for subsystem sizes of 2, 4, 6, respectively.System size is N = 10 and the coupling parameters are Ja = Je, Jx = 0.1Ja.

Figure
Figure S6.(a) Schematic of the action of the Hamiltonian on neighboring dimers.(b)Example of a sequence of states that is obtained using only allowed processes.The color of the singlet and triplet states indicates if they have been "exchanged" an even (blue) or odd (red) number of times.This exactly corresponds to the parity of the sum of the number of holons and doublons to their left.If one gets +1 for all blue triplets and red singlets, and −1 for all blue singles and red triplets, then for all states in the sequence the sum is the same and gives Q = −2, showing the conservation of this charge.

20 WDFigure S7 .
Figure S7.Distribution of energy level spacings for the model in Eq. (1) with N = 20 sites.Data is for half filling and Q = 0, with 25 disorder realizations.The level statistics displays excellent agreement with Wigner-Dyson statistics.The inset shows the average consecutive spacing ratio ⟨r⟩ for each realization, which is seen to be very close to 0.53, as expected for a chaotic system.Data is for Ja = 3 and J e,k ∈ [2, 2.5], ω k ∈ [0.5, 1.5] drawn from a uniform distribution.

Figure
Figure S8.Overlap between the |Π⟩, |ϕJ ⟩ and |ϕL⟩ states and the eigenstates of the model in Eq. (1).|Π⟩ only has overlap on the scarred eigenstates of the first family, while |ϕJ ⟩ and |ϕL⟩ have no overlap on them.For both of these states the scarred states of the second family dominate.Data is for N = 18, Ja = 4 and J e,k ∈ [4, 4.5], ω k ∈ [0.5, 1.5] drawn from a uniform distribution.

Figure S9 .
Figure S9.Fidelity and bipartite entanglement entropy over time following the quench from various initial states indicated in the legend.The thin black lines are randomly chosen Fock basis states at half filling.For |Π⟩ and |ϕJ ⟩ there is no visible growth of entanglement entropy, while for |ϕL⟩ the growth is strongly suppressed.Data is for N = 18, Ja = 4 and J e,k ∈ [4, 4.5], ω k ∈ [0.5, 1.5] drawn from a uniform distribution.1.0| 2 |Π |φ J |φ L

5 Figure
Figure S11.(a) Entanglement entropy computed from the analytical formula of the entanglement spectrum.The gray and black lines represent the expected large-M limit and they are in good agreement with the data already for N ≈ 20.(b) Difference between the data and the large-M prediction.For all three cases, the difference approximately decays as ∝ 1/M .

Figure
FigureS12.(a) Entanglement entropy of the scarred state of the second family with E = 0 for random realizations, for the analytical formulas, and from numerical minimization/maximization. As N increases, the analytical and numerical results converge.(b) Entanglement entropy of eigenstates in the Q = 0 sector for N = 12 in the parameter regime found to maximize the entanglement entropy of the scarred state of the second family.The non-scarred eigenstates concentrate around an arc, as is typical from chaotic systems.This shows that high entanglement entropy can be obtained without getting close to a fine-tuned integrable point where all parameters are identical.

40 Figure S13 .
Figure S13.Entanglement entropy of the eigenstates for a chain with N = 18, Ja = 4, J e,k ∈ [2, 6], and ω k ∈ [1, 3].Each subplot corresponds to a different sector of Q.The scarred states of the first and second families are denoted by pentagons and squares respectively.

k
Pd↔u involves a particle-hole exchange, the number of particles in subsystems 1 and 2 must sum to M .So by construction, rainbow scars in our model can only exist at half-filling.