Quantum dynamics of a fully-blockaded Rydberg atom ensemble

Classical simulation of quantum systems plays an important role in the study of many-body phenomena and in the benchmarking and verification of quantum technologies. Exact simulation is often limited to small systems because the dimension of the Hilbert space increases exponentially with the size of the system. For systems that possess a high degree of symmetry, however, classical simulation can reach much larger sizes. Here, we consider an ensemble of strongly interacting atoms with permutation symmetry, enabling the simulation of dynamics of hundreds of atoms at arbitrarily long evolution times. The system is realized by an ensemble of three-level atoms, where one of the levels corresponds to a highly excited Rydberg state. In the limit of all-to-all Rydberg blockade, the Hamiltonian is invariant under permutation of the atoms. Using techniques from representation theory, we construct a block-diagonal form of the Hamiltonian, where the size of the largest block increases only linearly with the system size. We apply this formalism to derive efficient pulse sequences to prepare arbitrary permutation-invariant quantum states. Moreover, we study the quantum dynamics following a quench, uncovering a parameter regime in which the system thermalizes slowly and exhibits pronounced revivals. Our results create new opportunities for the experimental and theoretical study of large interacting and nonintegrable quantum systems.


I. INTRODUCTION
Characterizing and understanding the emergent behavior of quantum many-body systems has been a key research focus of condensed-matter physics.Many questions still remain unanswered owing to the complexity of these systems.Recent experimental advancements have created new avenues to approach these open problems.It is now possible to experimentally study the dynamics of large quantum systems and to directly probe physical phenomena such as the growth of entanglement and the presence or absence of thermalization [1][2][3][4][5][6][7][8][9].The high level of precise control required for these experiments moreover constitutes a crucial step towards building large-scale quantum computers.
Ultracold atoms have been a particularly fruitful platform for this line of research [10][11][12][13][14].By means of cooling and trapping, one can prepare large ensembles of identical atoms.The internal degrees of freedom of each atom, such as distinct hyperfine levels, give rise to a discrete state space suitable for quantum information processing.While the states can be manipulated by external electromagnetic fields, they also retain coherence for a long time due to their weak coupling to the vacuum.The interaction between atoms in low-lying states is typically weak, which is beneficial for preserving the collective quantum state, but prevents the creation of complex states with large entanglement.
This limitation can be overcome, for example, by promoting atoms to highly-excited Rydberg states, which interact strongly via the van-der-Waals interaction [15].In this work, we investigate the quantum dynamics of The strong van der Waals interaction between two Rydberg states defines a blockade radius R b , within which the excitation of more than two atoms to the Rydberg state is effectively forbidden.We assume that R b is greater than the largest distance between any pair of atoms.Equivalently, a sphere of radius R b /2 around an atom in the Rydberg state (purple) is assumed to intersect the sphere of the same size centered at any other atom (orange).
the collective many-body system in the regime of the socalled Rydberg blockade, where the van-der-Waals interaction is so strong that it prevents excitation of more than one atom within a certain volume [16][17][18].In practice, a fully-blockaded ensemble may contain hundreds of atoms [19].When trapped individually in an array of optical tweezers, typical values of the blockade radius and the separation between atoms allow for a fully blockaded system comprising tens of atoms [20][21][22][23][24].The Rydberg blockade realizes an effective multi-atom interaction, which can be used to directly implement few-qubit gates without a decomposition into two-qubit gates [25][26][27][28][29][30].
To advance the development of quantum technologies arXiv:2311.18616v1[quant-ph] 30 Nov 2023 based on Rydberg atoms, it is desirable to develop classical simulations that enable efficient benchmarking.Owing to its nonlocality, the Rydberg blockade is expected to lead to rapid growth of entanglement.Many classical simulation techniques such as tensor networks will therefore be limited to short times or small systems.We address this challenge by considering a constrained model with a high degree of symmetry.In particular, we analyze instances where all atoms are within the blockade radius of each other, as shown by Fig. 1(b), and all driving fields are homogeneous.Due to the separation of scales inherent to the Rydberg blockade and thus insensitivity to the exact value of the interaction, the system is invariant under any permutation of the atoms.This symmetry enables numerical simulation of the quantum dynamics of ensembles containing several hundreds of atoms for arbitrarily long times.
Our formal approach is based on a representation theoretical analysis of the quantum Hamiltonian, which we introduce in section II.The analysis may be viewed as a generalization of Dicke states from two-level to threelevel atoms [31].Dicke states correspond to eigenstates of total angular momentum, which are mathematically constructed as the irreducible representations of the Lie group SU (2).For atoms with three levels, we have to instead consider the irreducible representations of SU (3).We apply this formalism in section III to construct pulse sequences that prepare arbitrary permutation-invariant quantum states.In section IV, we investigate the dynamics of a large system initialized in product states.We uncover a parameter regime in which the system exhibits robust revivals.The revivals indicate a slow progression to equilibrium despite the strong interactions.We show that the revivals are intimately linked to the symmetries of the model and provide a simple explanation in terms of an effective spin model.We summarize our findings in section V and discuss directions for future research.

A. System Hamiltonian
Throughout this work, we consider a set of n identical atoms whose dynamics are restricted to three internal states: two hyperfine states denoted by |0⟩ and |1⟩, and a highly excited Rydberg state |r⟩.Transitions between these states are controlled by two external driving fields as shown in Fig. 1(a).The first drive (Rabi frequency Ω 1 , detuning ∆ 1 ) acts on the |0⟩ ↔ |1⟩ transition, while the second drive (Ω 2 , ∆ 2 ) acts on |1⟩ ↔ |r⟩.Both drives are assumed to be global, meaning that the Rabi frequencies and detunings are the same for each atom.
The large polarizability of the Rydberg state gives rise to strong van der Waals interactions between nearby atoms, which creates an energy shift if two atoms are in the Rydberg state.The energy shift scales with the sixth power of the inverse distance [15].At short distances, the energy shift can therefore be much larger than any other scale, in which case states with more than one Rydberg excitation can be adiabatically eliminated.This mechanism is the so-called Rydberg blockade.The adiabatic elimination is valid at distances shorter than the blockade radius, R b , which depends on the driving strength [15].We assume that R b exceeds the largest separation between any pair of atoms (see Fig. 1(b)), so that it is a good approximation to assume that no more than one Rydberg excitation exists in the system.We will thus restrict our analysis to the subspace containing zero or one Rydberg excitation.The Hamiltonian describing the system can then be written, in the rotating frame and within the rotating-wave approximation, as where the subscript i indicates which atom an operator acts on and P projects onto the subspace with at most one Rydberg excitation.As a consequence of the global drive and the insensitivity of the Rydberg blockade to the distance, the Hamiltonian is invariant under permutations of the atoms.To take advantage of this symmetry, we introduce collective raising operators We further define the lowering operator T − = (T + ) † and the Hermitian operators T x = (T + + T − )/2, T y = (T + − T − )/(2i), and T z = [T + , T − ]/2 with analogous expressions for U and V .Together, these operators may be viewed as a generalization of collective spin operators to the setting of three-level systems.The Hamiltonian can be rewritten in terms of these collective operators as where we dropped an unimportant term proportional to the identity.

B. Review of irreducible representations of su(3)
The operators T α , U α , and V α (α ∈ {x, y, z}) generate the Lie algebra su(3) of the special unitary group SU (3).We can see this by considering a single atom, for which these operators are represented by traceless, Hermitian 3×3 matrices.Moreover, any traceless, Hermitian 3 × 3 matrix can be expressed as a linear combination of these operators with real coefficients.Hence, the algebra spanned by the operators corresponds to the defining representation of su(3).We note that the 9 operators are overcomplete since the Lie group SU(3) is only 8 dimensional.Indeed, one of the operators is redundant due to the linear relation The Lie algebra generated by a set of operators is fully determined by their commutation relations.We provide the full set of commutation relations between T α , U α , and V α in appendix B 2. Since the commutation relations are independent of the number of atoms, these operators also describe a representation of su(3) for more than one atom.The representation is reducible except for the special case of a single atom.Our strategy is to decompose the representation into a direct sum of irreducible representations, in which the operators T α , U α , and V α become block diagonal.We will see that the projector P , and hence the Hamiltonian, is block diagonal in the same basis.This greatly simplifies the computation of many quantities such as the eigenenergies.Formally, the block-diagonal form exists as a consequence of Schur-Weyl duality for permutation-invariant systems (see appendix A).In what follows, we pursue a less abstract viewpoint and focus on an explicit construction of the matrix representation of the irreducible blocks of H.
Before we turn to the irreducible representations of su(3), let us review the simpler case of su(2).This is particularly instructive because the operators T α generate a su(2) subalgebra of su(3).This follows from the fact that for a single atom, these operators act on the two-dimensional subspace spanned by {|0⟩ , |1⟩}.Alternatively, one can verify that they satisfy the familiar commutation relations of angular momentum, T α , T β = i γ∈{x,y,z} ε αβγ T γ , where ε αβγ is the Levi-Civita symbol.Similar arguments apply to the subalgebras generated by U α and V α .
The finite-dimensional, irreducible representations of su(2) are commonly labeled by the spin s, which is a non-negative integer multiple of 1/2.The dimension of the representation is 2s + 1.There exists an orthonormal basis {|m⟩}, where m runs from −s to s in integer increments.In this basis, the matrix elements of the generators of the Lie algebra, or spin operators, are fully determined by the relations T ± |m⟩ = s(s + 1) − m(m ± 1) |m ± 1⟩ and T z |m⟩ = m |m⟩.The basis states belonging to an irreducible representation are commonly illustrated using a so-called weight diagram [32,33], shown in Fig. 2(a).Each dot represents a basis state |m⟩, ordered from left to right by increasing m.Adjacent states are connected by nonvanishing matrix elements of T ± .
The above concepts set a stage for our discussion of the representations of su (3).We restrict ourselves to the aspects most pertinent to our problem and refer the The colored arrows show the directions in which the three distinct raising operators act.(c) In general, the weight diagram of the (p, q) representation of su(3) takes the form of a hexagon with side lengths p and q.Each solid dot and each surrounding circle correspond to a basis state.The basis states are simultaneous eigenstates of T z , U z , and V z .The respective eigenvalues mt, mu, and mv can be read off at the orthogonal projection onto the corresponding coordinate axis.The origin of the axis lies at the center of the diagram.(d) The matrix elements of the raising operators T + , U + , and V + for (p, q) = (0, 2) and the same color coding as in (b).
reader to references [32,33] for a more comprehensive exposition.To describe the irreducible representations, one can generalize the weight diagrams introduced above.For su(3), the nodes of the weight diagrams occupy a two-dimensional triangular lattice, in contrast to the onedimensional lattice in the case of su(2).The raising and lowering operators T ± , U ± , and V ± induce transitions between neighboring nodes along the three distinct lattice directions.The representation corresponding to the three basis states of a single atom hence form a triangle as shown in Fig. 2(b).
Two further examples of weight diagrams of irreducible representations of su(3) are shown in Fig. 2(c) and (d).In general, each irreducible representation forms a hexagon that can be labeled by two nonnegative integers (p, q).Three of the sides of the hexagon have length p, the other three length q.If either p = 0 or q = 0, the hexagon simplifies to a triangle (e.g., Fig. 2(b) and (d)).When p = q = 0, the weight diagram consists of a single node.The basis states associated with each node are indicated by black dots and circles.A single black dot means that there is one state associated with this node.Every circle surrounding a dot corresponds to an additional, degenerate basis state.The outermost nodes are always singly occupied.Moving inwards, the number of states per node increases by one if the nodes in the previous step formed a hexagon; if they form a triangle, the number of states stays the same.Following these rules, the (p, q) representation can be shown to have dimension (p + 1)(q + 1)(p + q + 2)/2 [32].
The basis states are simultaneous eigenstates of T z , U z , and V z .They are arranged in the weight diagram in such a way that the operators T ± , U ± , and V ± induce transitions between states occupying adjacent nodes along the three distinct lattice directions.Consequently, we can introduce three coordinate axes that determine the eigenvalues m t , m u , and m v respectively belonging to T z , U z , and V z at each point of the weight diagram as shown in Fig. 2(c).We can hence obtain the matrix elements of the diagonal operators T z , U z , and V z simply from the layout of the basis states in the weight diagram.
To compute the matrix elements of the raising and lowering operators T ± , U ± , and V ± , we make use of the fact that the states along a single line in a lattice direction form the basis of a representation of su (2).It follows that the matrix elements of the operators acting along these lines obey the rules of standard spin operators.However, care must be taken when constructing the matrix elements along different directions, especially when there is more than one state per node, to ensure that all commutation relations of su(3) are satisfied.We present a general recipe to compute the matrix elements relevant for our purposes in appendix B 3. As an illustration, Fig. 2(d) shows the matrix elements of the raising operators T + , U + , and V + in the (p, q) = (0, 2) representation.While the magnitude of the matrix elements can be understood by viewing each line along a lattice direction as either a spin-1/2 or spin-1 system, some matrix elements are inevitably negative because the operators T + , U + , and V + do not mutually commute.

C. Representation of n atoms
Following the above steps, we can construct explicit matrix representations of the operators T α , U α , and V α for any irreducible representation (p, q).In this section, we describe how the reducible representation of n atoms, defined by Eqs. ( 2)-( 4), decomposes into these irreducible representations.While p and q can in principle be any nonnegative integers, the values that arise for n atoms are restricted.This is similar to the addition of spins, where n spin-1/2 systems can combine into a total spin with quantum number of at most n/2.Given n atoms with three levels, the allowed values of p and q are in oneto-one correspondence with partitions of n into a sum of three non-negative integers [32,33].More concretely, let λ 1 ≥ λ 2 ≥ λ 3 ≥ 0 be three integers such that Each such partition of n corresponds to an allowed irreducible representation with In what follows, we will refer to an irreducible representation interchangeably by the labels (p, q) or by a partition λ = (λ 1 , λ 2 , λ 3 ).The Rydberg blockade places further constraints on the relevant irreducible representations by limiting the allowed basis states within each sector to have at most one Rydberg excitation.The number of atoms in the |0⟩, |1⟩, and |r⟩ states is well defined for every node in a weight diagram.This follows from expressing the number operators as together with the fact that the basis states are eigenstates of T z , U z , and V z .We note that the occupation numbers change along axes that are perpendicular to the m t , m u , and m v axes in Fig. 2(c).
The Rydberg blockade restricts the relevant Hilbert space to n r ≤ 1.We can therefore account for it by discarding all nodes in a weight diagram for which this condition is violated.Given the orientation of the m u and m v axes, Equation (10) implies that n r is constant along any horizontal row and increases by one for each row from top to bottom.Hence, the Rydberg blockade simply selects up to two horizontal rows from each diagram.Which rows are selected depends on the number of atoms n and on the irreducible representation (p, q).For the top-right corner of a weight diagram, we have U z = q/2 and V z = (p + q)/2, such that n r = (n − p − 2q)/3 = λ 3 in the top row.If λ 3 = 0, the top two rows are selected, which contain zero or one Rydberg excitation.If λ 3 = 1, then only the top row is kept.Irreducible representations with larger values of λ 3 contain only states with more than one Rydberg excitation and therefore do not participate in the dynamics.Below, we will often assume that the initial state contains no Rydberg excitations.In this case, it suffices to consider partitions with λ 3 = 0.
Let us illustrate these notions with the concrete examples of n = 2 and n = 3 atoms.For two atoms, there are two partitions: (λ 1 , λ 2 , λ 3 ) = (2, 0, 0) and (λ 1 , λ 2 , λ 3 ) = (1, 1, 0) corresponding to (p, q) = (2, 0) and (p, q) = (0, 1), respectively.The weight diagrams of these representations is shown in Fig. 3(a), where all states with more than one Rydberg excitation have been grayed out to reflect the blockade constraint.Explicit expressions for the Hamiltonians and the states in both the computational basis and the block-diagonal basis are given in appendix B 4. Similarly, Fig. 3(b) shows all possible partitions and the corresponding weight diagrams for n = 3 atoms.To account for the 27 basis states of the full Hilbert space (without blockade constraint), the partition (λ 1 , λ 2 , λ 3 ) = (2, 1, 0) occurs with multiplicity µ λ = 2.This is analogous to the fact that three spin-1/2 particles can form a system with total spin 1/2 in two different ways.
To further elucidate the role of the multiplicity of the irreducible representations, we write the quantum state of n atoms as where the sum of λ runs over the partitions of n.The multiplicity µ λ is accounted for by the index i, while the index j refers to the different states of the weight diagram.The number d λ is the dimension of the irreducible representation restricted by the Rydberg blockade.We provide explicit expressions for both µ λ and d λ in appendix B 1.
As a consequence of Schur-Weyl duality (see appendix A), the Hamiltonian, as well as any other operator that is invariant under permutation of the atoms, decomposes into blocks that depend on the irreducible representation but not on the multiplicity index i.Put differently, the Hamiltonian acts nontrivially only on the index j of the above basis.We may therefore express the Hamiltonian as where I µ λ is the µ λ × µ λ identity matrix.The matrix H λ can be constructed by inserting the λ representation of the operators T α , U α and V α into Eq.( 5) and restricting the Hilbert space according to Rydberg blockade constraint (see appendix B 3).Although we labeled the blocks H λ by the partition λ = (λ 1 , λ 2 , λ 3 ), they only depend on the parameters p and q of the irreducible representation and are independent of n.

D. Efficient computation
Let us briefly summarize the insights gained above.The Hamiltonian can be brought into a block-diagonal form, where each block is associated with an irreducible representation of SU(3) labeled by two integers p and q.The matrix elements of each block can be computed following the recipe in appendix B 3. Crucially, it is not necessary to compute the unitary transformation that brings the Hamiltonian into the block-diagonal form.This is an advantage from a computational perspective because the complexity of classically constructing this transformation increases exponentially with the number of atoms [34].However, the set of quantities that can be computed without the explicit basis transformation is restricted due to the limited knowledge about the relation between the basis states in the block-diagonal form and the original tensor-product basis.
The full energy spectrum of the Hamiltonian is an example of a physically relevant quantity that can be computed without knowledge of the basis states.We will show in Section IV that it is also possible to compute the number of atoms in the |0⟩, |1⟩, or |r⟩ state after evolving an initial product state for an arbitrary time t.In both cases, the computational effort scales polynomially with the number of atoms, n.To see this, we recall that the Rydberg blockade constraint requires that λ 3 ≤ 1, which implies that the number of distinct irreducible representations increases linearly with n.The dimension d λ of the unblockaded subspace in a given irreducible representation increases linearly with p (see appendix B 1), which is in turn bounded from above by n.Since the eigenvalues and eigenvectors of a matrix of size m × m can be determined in time O(m 3 ) [35], the computational effort to exactly diagonalize all irreducible representations is O(n 4 ).
The polynomial scaling represents an exponential improvement over brute-force diagonalization of the Hamiltonian in the blockaded Hilbert space, which, considering the Rydberg blockade, has dimension 2 n + n2 n−1 .We remark that the exponential size of the Hilbert space is accounted for by large multiplicities µ λ (see appendix B 1).
The computational speedup is enabled by the fact that it suffices to diagonalize each irreducible representation once, regardless of its multiplicity.
We note that the abstract formalism in terms of irreducible representations can be useful even if knowledge of the basis states is required.We discuss such a situation in the next section, where we explore the preparation of quantum states that are invariant under permutations.Efficient computation is possible for this special case because all permutation-invariant states of n atoms are part of the (n, 0) representation.

III. PREPARATION OF PERMUTATION-INVARIANT QUANTUM STATES
Having established our formalism to describe a fully blockaded system of n atoms, we apply it to the preparation of arbitrary permutation-invariant quantum states.First, we note that all permutation-invariant quantum states belong to the (p, q) = (n, 0) representation, which includes the states where all atoms are in the |0⟩ or |1⟩ state.This follows formally from the Schur-Weyl duality (see appendix A).The basis states in the top row of the (n, 0) representation are the so-called Dicke states |D n n0 ⟩ [31], which consist of equal superposition of n 0 atoms in the |0⟩ state and n 1 = n − n 0 atoms in the |1⟩ state, where n 0 ranges from 0 on the left to n on the right of the top row.Any permutation-invariant state in the |0⟩ and |1⟩ basis can be decomposed into a superposition of Dicke states.We show below that one can design a pulse sequence to prepare an arbitrary permutation-invariant state starting with all atoms in the |0⟩ states.
We explain the protocol with a simple example of preparing a two-qubit Greenberger-Horne-Zeilinger (GHZ) state (or Bell state), (|00⟩ + |11⟩)/ √ 2, as illustrated in Fig. 4. The protocol closely follows the steps originally proposed in Ref. [17].The key idea is to derive a pulse sequence for the reverse process of preparing the |00⟩ state from the GHZ state.The actual state preparation can then be obtained by time reversal of this pulse sequence.
Starting from the GHZ state, one can progressively move the population in the top-left node of the weight diagram (|11⟩ state) to the top-right node (|00⟩ state), as shown in Fig. 4. In each of the steps, we choose a pulse such that the population of the left-most occupied node is moved to the right.To this end, we alternate between pulses that act on the |1⟩ ↔ |r⟩ and on the |0⟩ ↔ |r⟩ transition.The pulse on the |1⟩ ↔ |r⟩ transition can be realized using then Rabi drive Ω 2 in the Hamiltonian, Eq. (5).By contrast, we do not allow for a direct drive of the |0⟩ ↔ |r⟩ transition in the Hamiltonian.It is possible, however, to effectively realize such a pulse by first applying a π-pulse on the |0⟩ ↔ |1⟩ transition, which exchanges the |0⟩ and |1⟩ state of each atom.The desired pulse can then be applied to the |1⟩ ↔ |r⟩ transition, be- Owing to the Rydberg blockade, the pulses connect pairs of states between the first and second row of the weight diagram.Hence, we can determine the pulse that empties the population of the left-most state by exactly solving the dynamics of a two-level system.As we can see from Fig. 4(a) and (b), the first two steps are π-pulses, transferring the population from the top-left node to the bottom-left node and then to the top-middle node.A πpulse can be implemented by setting all detunings to zero and choosing the pulse duration t and the appropriate Rabi frequency Ω such that kΩt = π.The numerical factor k accounts for the matrix element of U ± or V ± connecting the two states.Following the discussion of section II B, it is given by k = √ 2 for the first pulse and k = 1 for the second pulse.
The state-dependent matrix elements lead to a partial transfer of the population from the top-right node (|00⟩ state) to the bottom-right node during the second pulse.Therefore, the third pulse is not a π-pulse.Nevertheless, we can readily determine the required pulse numerically by keeping track of matrix elements and the amplitude of each basis state, as illustrated in Fig. 4. The approach readily generalizes to more than 2 atoms.The numerically calculated pulse sequences for preparing GHZ states for n = 2, 3, and 4 atoms are included in appendix C.
We note that the existence of the Rydberg energy level |r⟩ and the blockade effect are crucial for this protocol by allowing us to isolate effective two-level systems.If we were to directly drive the |0⟩ ↔ |1⟩ transition, all states on the top row would evolve together, rendering it impossible to completely move the population from the GHZ state to the state where all atoms are in the |0⟩ state.This is expected as one cannot create an entangled state from a product state using only the noninteracting hyperfine states.
The above protocol straightfowardly generalizes beyond GHZ states to arbitrary permutation-invariant states, including the W state.In general, it consists of n pulses on the |1⟩ ↔ |r⟩ transition alternating with n pulses on the |0⟩ ↔ |r⟩ transition.The latter are in practice each decomposed into 3 pulses, resulting in 4n pulses in total.We highlight that the linear dependence of the number of pulses on the system size is optimal [36].The parameters of each pulse can be computed numerically considering only effective two-level systems such that the computational cost increases only linearly with the system size.

A. Method
We now use our formalism to study the equilibration of a permutation-invariant Rydberg system.We emphasize that we only consider the unitary dynamics of a closed system.In contrast to an open system with dissipation [37], the system remains in a pure state at all times and does not converge to a mixed steady state.Equilibration in a closed system occurs as the off-diagonal elements of an observable in the energy basis dephase.At late times, the time-averaged expectation value of the observable converges to the average of its diagonal elements in the energy basis, weighted by the overlap of the initial state with the corresponding eigenstate [38].
We numerically explore the time evolution of an observable for different initial states.In the examples below, we focus on the number of Rydberg excitations, n r , although the approach readily applies to any observable that is invariant under permutations of the atoms.As pointed out in Section II C, any such observable A can be written as where the direct sum runs over the irreducible representations of SU(3) labeled by λ.The matrices A λ describe the action of A on the vector space associated with the irreducible representation.In the case of n r , these irreducible blocks can be readily constructed using Eq.(10).Given a state |ψ⟩, the expectation value of A may be expressed as where ρ λ is the (unnormalized) state obtained by projecting |ψ⟩ ⟨ψ| onto the irreducible representation labeled by λ and tracing out its multiplicity.Using the notation of Eq. ( 11), ρ λ can be written explicitly as Here, the states {|λ, j⟩} form an orthonormal basis of the d λ -dimensional vector space associated with the irreducible representation.Since the projection commutes with the Hamiltonian, we can apply the time evolution to each ρ λ separately.Given the initial state ρ λ (0), we can classically compute the time-evolved state ρ λ (t) = e −iH λ t ρ λ (0)e iH λ t at a computational cost that increases polynomially with the system size.Hence, it is possible to efficiently determine the time-dependent expectation value of permutationally invariant observables.We focus on initial states with a well-defined number of atoms in the |0⟩ and |1⟩ states and no atom in the Rydberg state, |r⟩.According to Eqs. ( 8)-( 10), such states are eigenstates of T z , U z , and V z .The respective eigenvalues uniquely specify the node in the weight diagram of each irreducible representation.Moreover, since n r = 0, only irreducible representations with λ 3 = 0 contribute and the node lies in the top row.Because all nodes in the top row host nondegenerate states, the occupation numbers specify a unique initial state for each irreducible representation.We denote this state by |λ, n 0 , n 1 ⟩, where n 0 and n 1 are the number of atoms in the |0⟩ and |1⟩ states, respectively.Taking this state to be normalized, we have where w λ is the probability that the system occupies the irreducible representation labeled by λ.We show in appendix B 5 that if λ 3 = 0 and |n 0 − n 1 | ≤ p ≤ n, and w λ = 0 otherwise.The representations A λ with A = n r are diagonal matrices with entries 0 for the states corresponding to the top row of the weight diagram and 1 for the ones in the bottom row.

B. Numerical results
We next numerically study the dynamics of large systems with different initial states.Figure 5 shows the results for n = 100 atoms with initial product states of the form |ψ(0)⟩ = |0⟩ ⊗n0 |1⟩ ⊗n1 .We remark that due to permutation invariance, it is irrelevant which of the atoms are in the |0⟩ state and which in the |1⟩ state.The system is subject to simultaneous driving of the |0⟩ ↔ |1⟩ and the |1⟩ ↔ |r⟩ transitions with equal Rabi frequencies, Ω 1 = Ω 2 = Ω.The field driving the hyperfine transition is detuned by ∆ 1 = Ω, whereas the Rydberg drive is resonant, ∆ 2 = 0. We note that only the (p, q) = (n, 0) representation is relevant for Figs.5(a) and (c), where all atoms occupy the same state, rendering the state invariant under permutations.This is not the case in Fig. 5(b), for which we have to sum over all irreducible representations satisfying λ 3 = 0.For all three initial states in Fig. 5, the Rydberg population oscillates rapidly.We note that the typical period of the oscillation is much shorter than 2π/Ω 2 , which is the duration of a 2π pulse between the |1⟩ and |r⟩ states of a single atom.This can be explained by the enhanced Rabi frequency experienced by states with more than one atom in the |1⟩ state.The Rabi frequency experienced by the |1⟩ ⊗n state is given by √ nΩ 2 , which shortens the period by a factor 1/ √ n.Other states experience a weaker enhancement, which depends on both the number of atoms in the |0⟩ and |1⟩ states as well as the irreducible representation.Following the construction of matrix elements in Section II and appendix B 3, we find that the largest enhancement factor in the (p, q) representation is (n + p)/2.Hence an enhancement factor proportional to √ n is expected independent of the initial state.In addition to the fast oscillation, the dynamics exhibit a much slower modulation of the oscillation amplitude.When the atoms all start in the state |0⟩ (Fig. 5(a)), the fast oscillations as well as the slow modulation possess little structure.By contrast, when half or all atoms start in the |1⟩ state (Figs.5(b) and (c)), the system undergoes damped oscillations with strong revivals after a time of more than 100/Ω.Revivals in such a large system are surprising and are often indicators of symmetries such as the permutational invariance in the present setting.
We characterize the revivals by determining the time and strength of the revival.We accomplish this numer-ically by sorting the evolution times according to their corresponding Rydberg population.The largest Rydberg populations occur at short times.As the Rydberg population decreases, the sorted evolution time jumps up when the highest point of a revival is reached.Because the strength of the revivals is monotonically decreasing, this jump reliably identifies the first revival.We refer to the time and Rydberg population at this point as the revival time t rev and the revival strength ⟨n r ⟩ rev , respectively (see Fig. 5(b)).We plot both quantities in Figure 6 as a function of the number of atoms for the two initial states that exhibit pronounced revivals.
Figure 6(a) shows that the time of the first revival increases with √ n.The functional dependence of the revival strength on n is not evident from the data in Fig. 6(b), but the revival strength remains above 0.9 even for n = 300 atoms.This indicates that the first revival stays pronounced in large systems even as it shifts to later times.The strength of subsequent revivals, however, decreases and data at very late times show that they eventually disappear (see appendix D 2).

C. Spin model
To shed light on the physical origin of the revivals, we consider the dynamics of the initial state |1⟩ ⊗n in more detail.In this case, the evolution is restricted to the irreducible representation (p, q) = (n, 0), whose weight diagram is shown in Fig. 7(a).The top row of this diagram, where no atom is in the Rydberg state, may be viewed as the states of a spin with total angular momentum quantum number s = n/2.Similarly, the bottom row, where exactly one atom is in the Rydberg state, represents a spin s = (n − 1)/2.In the limit of large n, we may ignore the difference between these spins.We formalize this idea by adding an unphysical state to the bottom row.This additional state will have little impact on the dynamics provided its amplitude remains negligible throughout the evolution.We expect this to be indeed the case because the initial state is situated at the top-left corner of the weight diagram and the detuning ∆ 1 = Ω creates an energy difference of order n∆ 1 between opposite ends of the diagram such that occupying the added state is energetically unfavorable.The addition of the unphysical state leads to a great simplification because both rows now have n + 1 states.This allows us to factorize the Hilbert space into a spin s = n/2 and a two-level system whose states |0⟩ and |1⟩ correspond to the number of atoms in the Rydberg state.As a further approximation, we modify the matrix elements of T x,y,z in the bottom row, which originally correspond to those of a spin s = (n − 1)/2.In the limit of large n, it is justified to replace the matrix elements by those of a spin s = n/2.The spin-model Hamiltonian resulting from these approximations can be compactly expressed as where σ x,y,z are Pauli matrices and S x,y,z are the standard spin operators with total angular momentum s = n/2.For simplicity, we assumed that Ω 1 and Ω 2 are real.The operators S x,y,z couple the same states as T x,y,z with the above modification of the matrix elements and the additional, unphysical state.The term sI 2s+1 − S z accounts for the matrix elements of the U ± operators in the original Hamiltonian (see Fig. 7(a)).We note that the Pauli matrices are related to the occupation of the Rydberg state by n r = I 2s+1 ⊗ (I 2 − σ z )/2.In the special case ∆ 2 = 0, the states |±⟩ = (|0⟩ ± |1⟩)/ √ 2 of the two-level system are decoupled.It is therefore convenient to express the state of the system as The states |ϕ ± (t)⟩ evolve under the Hamiltonian In terms of this ansatz, the number of atoms in the Rydberg state is given by This implies that 1/2 The absolute value of the overlap between |ϕ + (t)⟩ and |ϕ − (t)⟩ hence determines the envelope of the Rydberg population, whereas the relative phase is responsible for the fast oscillations.
To verify the validity of the spin model, we plot in Fig. 7 ⊗n and the Hamiltonian parameters are the same as in Figs. 5 and 6.There is excellent agreement between the exact evolution and the envelope.We highlight that the spin model also captures the fast oscillations, which have been omitted from the figure for clarity.
The spin model enables us to explain the characteristics of the revivals displayed in Fig. 6.In fact, it is sufficient to work within a semiclassical approximation, which is valid in the limit of large s.We write S z = ⟨S z (t)⟩ + δS z (t) and expand sI 2s+1 − S z to linear order in δS z (t).This yields The expansion is justified if the fluctuation δS z (t) is small compared to s − ⟨S z (t)⟩.This is indeed the case for sufficiently large detuning ∆ 1 : Since ⟨S z (0)⟩ = −s, the detuning ensures that s − ⟨S z (t)⟩ is of order s at all times.
In comparison, the fluctuations in the initial state are smaller by a factor 1/ √ s according to ⟨δS z (0) 2 ⟩ ∼ √ s.We note that the uncertainty in the spin vector remains constant within the approximations of Eq. ( 22).Although higher-order corrections may cause the fluctuations to grow, we expect those terms to only become relevant at very late times.
By substituting Eq. ( 22) into the expression for the spin model Hamiltonian in Eq. ( 20), we find that the states |ϕ ± (t)⟩ evolve according to The term proportional to the identity leads to a phase difference between the states and thus causes oscillations in the Rydberg population with a frequency of order √ nΩ 2 .The second row of the Hamiltonian describes an effective magnetic field around which the states |ϕ ± (t)⟩ precess.Since there are no nonlinear terms, the states remain spin-coherent states at all times and can be simply described by the direction of their spin vector.The state overlap can be expressed as | ⟨ϕ + (t)|ϕ − (t)⟩ | = cos 2s (θ/2), where θ is the angle between the two respective spin vectors [39].The magnitude of the magnetic fields experienced by the two states differs by a term proportional to Ω 2 / s − ⟨S z (t)⟩.Hence, the two spin vectors separate and rephase over a time of order √ s, which explains the revival timescale t rev ∼ √ n observed in Fig. 6(a).The revival is imperfect because the rotation axes fluctuate such that the spins do not return to the exact initial direction but are displaced by an angle of order ∆θ ∼ 1/ √ s.We expect that the strength of the first revival, ⟨n r ⟩ rev , is independent of the system size since cos 2s (c/ √ s) tends to e −c 2 /2 in the limit of large s for any constant c.This is consistent with the weak dependence on the system size observed in Fig. 6(b).
We have shown that the spin model captures the salient features of the revivals when the system starts in the state |1⟩ ⊗n .The system fails to equilibrate quickly because it is well approximated by a linear spin model.As pointed out above, the revivals decay over time, which is also true within the spin model.The approximations of the spin model, however, render it inapplicable at very late times, as can be seen in Fig. 9 of appendix D 2.
The spin model can be adapted to the initial state |0⟩ ⊗n (see appendix D 1).The additional, unphysical state must be added to the left end of the weight diagram to ensure its occupation remains small.However, the resulting Hamiltonian is more complicated than Eq. ( 18) because driving the |1⟩ ↔ |r⟩ transition results in a term that simultaneously changes the large spin s and flips the two-level system.The dynamics can therefore not be separated as in Eq. (19).As a consequence, the large spin and the two-level system do not periodically disentangle, which explains the absence of revivals for this initial state.For the initial state |0⟩ ⊗n/2 |1⟩ ⊗n/2 , a similar analysis is more involved because it requires summing over multiple irreducible representations.Moreover, a two-level system coupled to a large spin is insufficient to capture the two states per node that appear in the second row of many representations.The pronounced revivals for this initial state are thus surprising and a more refined model is required to explain their physical origin.

V. SUMMARY AND OUTLOOK
In summary, we presented a general formalism to efficiently compute the quantum dynamics of an ensemble of three-level systems evolving under a permutationinvariant Hamiltonian.The high degree of symmetry allows us to numerically simulate ensembles of 100s of atoms that interact via the Rydberg blockade mechanism.We applied the formalism to derive sequences of at most 4n pulses to prepare arbitrary symmetric states of n atoms, including important entangled states such as the W state and the GHZ state.In addition, we studied the dynamics of an ensemble after a quench and uncovered a regime of parameters in which the system approaches equilibrium surprisingly slowly despite the strong nonlinearity due to the Rydberg blockade.We explain this slow equilibration and the appearance of revivals in terms of an effective spin model for permutation-invariant initial states.The physical origin of the revivals for other initial states will be the subject of future research.
Our formalism generalizes the notion of Dicke states from two-level to three-level systems.The approach is general and can accommodate other forms of interaction besides Rydberg blockade, provided the Hamiltonian remains invariant under permutations.It is possible to generalize our results beyond three-level systems to arbitrary d-level systems as the underlying Schur-Weyl duality also applies to the group SU(d).Future work may further extend the formalism to open systems governed by a permutationally invariant Lindblad operator [37,40].irreducible representation (p, q) by |r; t, m t ⟩.Here, r is a non-negative integer assigned to each row, with r = 0 corresponding to the top row.We further note that the states in the weight diagram are simultaneous eigenstates of U z and V z .Using the relations U z − V z + T z = 0 and r = (p + 2q − 2U z − 2V z ), we find It remains to determine the matrix elements of U ± and V ± .These matrix elements are more complicated because they couple states with different values of t.As we are only interested in matrix elements between states with at most one Rydberg excitation, we will restrict the remaining discussion to the top two rows of any weight diagram, which limits the range of values of t.The top row is always singly occupied with t = p/2.To determine the allowed values for t in the second row, we distinguish between four cases.p = q = 0: The irreducible representation is onedimensional and there is only a single state, for which t = 0. p = 0, q ̸ = 0: All nodes in the second row are singly occupied with t = 1/2.
The matrix elements can now be constructed from the commutation relations of the Lie algebra [32].Given an integer 0 ≤ k ≤ p/2, they can be concisely written as The matrix elements of U + and V + follow by taking the Hermitian conjugate.The expressions hold in all four of the above cases.When p = 0 or q = 0, some of the coefficients on the right-hand side vanish such that only states with either t = (p + 1)/2 or t = (p − 1)/2 are present.

Hamiltonian for n = 2 atoms
In this appendix, we give explicit expressions for the Hamiltonian for n = 2 atoms, both in the original basis and in the block-diagonal form.Taking into account the Rydberg blockade, the Hilbert space is spanned by the eight states {|00⟩ , |01⟩ , |0r⟩ , |10⟩ , |11⟩ , |1r⟩ , |r0⟩ , |r1⟩}.Using this ordering of the states, the Hamiltonian in Eq. ( 1) has the matrix representation Following the discussion of section II, this Hamiltonian can be brought into block-diagonal form with two blocks corresponding to the partitions λ = (2, 0, 0) and λ = (1, 1, 0).The blocks are explicitly given by Here, we labeled the blocks by (p, q) instead of the partition λ to highlight that their form also applies to other values of n.
One can explicitly check that H (2,0) and H (1,0) correspond to the matrix elements of H in the subspaces respectively spanned by These are the basis states associated with the nodes of the weight diagrams labeled as in appendix B 3. We observe that the states in the (p, q) = (2, 0) irreducible representation are invariant under permutation.The states in the (p, q) = (0, 1) representations change sign under permutation.The states can be constructed by first constructing the eigenstates of total angular momentum of two atoms with two levels |0⟩ and |1⟩.The remaining states can subsequently be obtained by applying the operators U − and V − .

Projector
In this appendix, we derive a closed-form expression for the overlap of a given state with the (p, q) representation.We consider initial product states of the form |s 1 ⟩ |s 2 ⟩ • • • |s n ⟩, where s i ∈ {0, 1} for all i ∈ {1, 2, . . ., n}.Since there are no Rybderg excitations for this initial state, the overlap with the (p, q) representation is equal to the probability that the product state carries total spin quantum number t = p/2 with respect to the operator T 2 .
Given a finite-dimensional representation ρ of a compact group G, the projector, P λ , onto an irreducible representation ρ λ of dimension d λ can be written as (see chapter 5 of reference [42]) Here, the integral runs over the Haar measure on G and χ λ (U ) = tr[ρ λ (U )] is the character of the group element in the irreducible representation.Since the Rydberg level plays no role, the group of interest is SU (2).The elements of SU(2) can be described by the rotation axis n = (sin θ cos φ, sin θ sin φ, cos θ) and the rotation angle α.The initial representation is defined by applying the same rotation to each atom, (B31) The dimension of the irreducible representation with total angular momentum t is d t = 2t + 1 and the character is given by (see example 12.23 of reference [33]) sin [(2t + 1)α/2] sin(α/2) .
We next apply these expressions to a product state with n 0 atoms in the state |0⟩ and the remaining n 1 atoms in the |1⟩ state.This yields ⟨ψ|P t |ψ⟩ = 2t + 1 2π  The integral over cos θ can be carried out by parts, resulting in a sum of elementary integrals over α.Assuming that n 0 ≥ n 1 , we arrive at ⟨ψ|P t |ψ⟩ = (2t + 1) The sum is zero unless 0 ≤ n 2 − t ≤ n 1 or, equivalently, π flip in the phase of the drive, which reduces the rotation angle to π − θ.

FIG. 1 .
FIG. 1. System setup.(a) We consider atoms with two hyperfine ground states |0⟩ and |1⟩ and a highly excited Rydberg state |r⟩.The hyperfine states are coupled by a driving field with Rabi frequency Ω1 and detuning ∆1.The state |1⟩ is further coupled to |r⟩ with Rabi frequency Ω2 and detuning ∆2.(b)The strong van der Waals interaction between two Rydberg states defines a blockade radius R b , within which the excitation of more than two atoms to the Rydberg state is effectively forbidden.We assume that R b is greater than the largest distance between any pair of atoms.Equivalently, a sphere of radius R b /2 around an atom in the Rydberg state (purple) is assumed to intersect the sphere of the same size centered at any other atom (orange).

FIG. 2 .
FIG. 2. Irreducible representations of su(2) and su(3).(a) Weight diagram of the spin s = 3/2 representation of su(2).Each dot corresponds to a basis state |m⟩.The raising operator T + induces transitions between the states as indicated by the orange arrows.The numbers next to the arrows give the associated matrix element.(b) The three levels of a single atom form the (p, q) = (1, 0) representation of su(3).The colored arrows show the directions in which the three distinct raising operators act.(c) In general, the weight diagram of the (p, q) representation of su(3) takes the form of a hexagon with side lengths p and q.Each solid dot and each surrounding circle correspond to a basis state.The basis states are simultaneous eigenstates of T z , U z , and V z .The respective eigenvalues mt, mu, and mv can be read off at the orthogonal projection onto the corresponding coordinate axis.The origin of the axis lies at the center of the diagram.(d) The matrix elements of the raising operators T + , U + , and V + for (p, q) = (0, 2) and the same color coding as in (b).
FIG. 3. (a) Irreducible representations arising for n = 2 atoms.The axis shows the number of Rydberg atoms.States with more than one Rydberg excitation are grayed out as they violate the blockade constraint.The table shows the partition (λ1, λ2, λ3), the labels (p, q), and the multiplicity µ λ associated with each weight diagram.(b) The same as (a) for n = 3 atoms.

FIG. 4 .
FIG.4.Preparation of the two-atom GHZ state starting from the initial state |00⟩.The preparation corresponds to the time reversal of the five pulses (a)-(e).The numbers next to the arrows indicate the product kΩt for each pulse, where Ω is the Rabi frequency for a single atom, t the duration of the pulse, and k the magnitude of the matrix element of the collective transition.The amplitude of each basis state (up to a global phase) is indicated by the numbers next to the nodes, which are shaded according to the corresponding occupation.

FIG. 5 .
FIG.5.Expected number of Rydberg excitations in a fully-blockaded ensemble of n = 100 atoms.The dynamics are shown for three different initial states, indicated at the top of each panel.The system is simultaneously driven by two external fields with Rabi frequencies Ω1 = Ω2 = Ω and detunings ∆1 = Ω and ∆2 = 0 (see Fig.1).The plots in the top row show the rapid oscillations of the Rydberg population for times between Ωt = 0 and Ωt = 20.

FIG. 6 . 5 .
FIG. 6.(a) Revival time trev and (b) revival strength ⟨nr⟩rev as a function of the number of atoms, n.Both quantities are determined by locating the largest value of ⟨nr⟩ at the first revival.The Hamiltonian parameters are the same as in Fig. 5.The solid curves in (a) are fits to the function a √ n+b.

FIG. 7 .
FIG. 7. (a) Illustration of the spin model approximation to the (n, 0) representation.By adding an unphysical state (gray circle connected to dashed line), the Hilbert space can be viewed as the tensor product of a spin s = n/2 and a two-level system σ.The colored lines indicate the transitions that are being driven.The purple numbers show the matrix elements for the |1⟩ ↔ |r⟩ transition.(b) Revivals of an ensemble of n = 300 atoms initialized in the state |1⟩ ⊗n .The Hamiltonian parameters are the same as in Fig. 5.The gray curve shows the exact values of the expectation value of the number of Rydberg excitations.The blue envelope was obtained from the approximate spin model.
The value of r increases by one each time that we move down by one row.The other two labels are angular momentum eigenvalues such thatT 2 |r; t, m t ⟩ = t(t + 1) |r; t, m t ⟩ , (B13) T z |r; t, m t ⟩ = m t |r; t, m t ⟩ , (B14)whereT 2 = (T x ) 2 + (T y ) 2 + (T z ) 2 .We may view m t as a horizontal coordinate in the weight diagram.Its values are restricted to {−t, −t + 1, . . ., t}.The total angular momentum quantum number t serves to resolve different states on multiply occupied nodes of the weight diagram.Since the operators T α form a su(2) subalgebra, we further have T ± |r; t, m t ⟩ = t(t + 1) − m t (m t ± 1) |r; t, m t ± 1⟩ .(B15)

TABLE I .
Pulse parameters to prepare a GHZ state of n = 2, 3, and 4 atoms.The detuning is zero for each pulse and the phase of drive is constant.The pulses are ordered to prepare the GHZ state from the |0⟩ ⊗n atoms, which is opposite to the order in which the parameters were computed.