Thermalization of a strongly interacting 1D Rydberg lattice gas

When Rydberg states are excited in a dense atomic gas the mean number of excited atoms reaches a stationary value after an initial transient period. We shed light on the origin of this steady state that emerges from a purely coherent evolution of a closed system. To this end we consider a one-dimensional ring lattice, and employ the perfect blockade model, i.e. the simultaneous excitation of Rydberg atoms occupying neighboring sites is forbidden. We derive an equation of motion which governs the system's evolution in excitation number space. This equation possesses a steady state which is strongly localized. Our findings show that this state is to a good accuracy given by the density matrix of the microcanonical ensemble where the corresponding microstates are the zero energy eigenstates of the interaction Hamiltonian. We analyze the statistics of the Rydberg atom number count providing expressions for the number of excited Rydberg atoms and the Mandel Q-parameter in equilibrium.

Rydberg atom is excited. One observes an initial (transient) phase with large contrast oscillations followed by a steady state. These features do not depend on the actual details of the system, i.e. the exact atomic interaction, the atomic arrangement and the boundary conditions. The data here is shown for a ring lattice with 25 sites. Each site is occupied by the same number of atoms N 0 and the atoms are resonantly coupled to a Rydberg state by a laser with Rabi frequency Ω 0 . The interaction strength between neighboring Rydberg atoms is taken to be infinite. The next-nearest neighbor interaction is zero. The time is given in units of the inverse collective Rabi frequency Ω −1 .

Introduction
In Rydberg states of alkali atoms the valence electron occupies a high lying orbit and is only loosely bound to the ionic core [1]. This large displacement between the charges (∼ 100 nm) gives rise to a large atomic polarizability [2] and strong interatomic interactions. This has been impressively demonstrated in an ultracold gas where atoms have been laser excited to Rydberg states from their electronic ground state [3]. The strong state-dependent interaction forbids the simultaneous excitation of two atoms when they are too close. This effect -the Rydberg blockade -has potential application in quantum information processing [4][5][6][7][8][9], the preparation of many-particle states [10] in optical lattices and the creation of non-classical light sources [11,12].
A coherently laser-driven gas of Rydberg atoms constitutes an ideal specimen to study the complex many-body dynamics of a strongly interacting quantum system. The underlying Hamiltonian can be formulated as a spin model [10,[13][14][15][16][17] whose properties are tunable through the laser parameters and the choice of particular Rydberg states. The experimental realization of such a system has been achieved in several labs and the detection of Rydberg atoms can be carried out with high efficiency. Hence, observables such as the number of Rydberg atoms can be recorded with high accuracy. A typical experiment starts with all atoms being initialized in the ground state. The laser which couples the ground to the Rydberg state is turned on for a certain time t after which the number of Rydberg atoms is detected [3,18,19]. The resulting signal is determined by the interplay between the laser driving and the strong state-dependent interaction [3].
The typical appearance of the temporal evolution of the Rydberg number n r obtained from theoretical studies [10,13,15,16] is depicted in fig. 1. At small times (t ≤ t trans ) n r rises first quadratically in t showing subsequently a number of oscillations with high contrast. After this transient period the oscillations diminish and n r approaches a constant value around which it performs small amplitude fluctuations which decrease with increasing system size. This general appearance is independent of the exact details of the system, e.g. the actual arrangement of the atoms and the boundary conditions. In experiments which are carried out in disordered gases the individual oscillations might not be observable since averaging over many experimental realizations also means averaging over many different spatial distributions of the atoms [14]. Each of these distributions gives rise to a different shape of the oscillating features in fig. 1 and eventually the oscillations are washed out.
The steady state shown in fig. 1 is the result of a purely coherent dynamics of a closed system and does not come about due to dissipation stemming from the coupling to an external bath. In experiments a fundamental cause of dissipation would be radiative decay of atoms which typically occurs on a timescale τ rad ∼ 10 − 100 µs. This is long compared to the time it takes to carry out the experiment τ exp ∼ 1µs. Consequently, it is well-justified to neglect radiative decay in theoretical models of gases of Rydberg-excited atoms when the time-evolution is restricted to time intervals smaller than τ rad . Since any numerical treatment can only encompass a finite number of degrees of freedom, there is always a revival time t rev after which the system returns to its initial state. The steady state which is observed in fig. 1 can thus be itself only a transient effect. However, in practice we have τ rad ≪ t rev and thus the common theoretical models which consider only coherent dynamics are invalid for such long times. Thus, the steady state we are referring to in this work exists in the time interval t trans < t < τ rad .
The purpose of this work is to study the origin of this steady state which occurs as result of a purely coherent dynamics of a closed system. This subject is closely related to the discussion about how and why a closed system which is prepared in a pure state actually thermalizes. I.e. why and how does such system assume a state in which the mean values of macroscopic observables are stationary and why can these mean values be calculated from the microcanonical ensemble? These questions are of fundamental interest [20,21] and have been investigated in a number of systems [22][23][24][25][26].
In the present work it is particularly insightful to study the evolution of the atomic gas in excitation number space, i.e. between subspaces of the system's Hilbert space which contain the same number of Rydberg atoms. Throughout, we employ the perfect blockade model of a Rydberg gas which has been extensively discussed in Refs. [13,16]: Two-level atoms (ground state |g , Rydberg state |r ) are confined to a ring lattice with N sites. Each site contains the same number (N 0 ) of atoms, which occupy a single spatial mode and no tunneling between adjacent sites occurs. I.e. the external degrees of freedom are frozen out. A laser with Rabi frequency Ω couples |r and |g , resonantly. The interaction between Rydberg atoms is so strong that the interaction energy of two atoms occupying the same site or adjacent sites is infinite ‡. The interaction between Rydberg atoms being separated by two or more sites is zero.
The physical degrees of freedom are the two internal states of the atoms located at individual sites, i.e. their measurement eventually determines n r . The state-dependent interaction in conjunction with the laser coupling, however, leads to such a strong mixing that perturbation theory in these physical degrees of freedom becomes meaningless. Instead, the eigenstates of the system are collective, delocalized excitations which contain no well-defined number of Rydberg atoms. When studying the evolution in excitation number space under the action of the laser, we find that the strong interaction causes quasi-random couplings between regions of the Hilbert space which contain a different (and well-defined) number of Rydberg atoms. This randomness allows us to derive an effective equation for the time-evolution of the probability of being in a subspace with certain fixed excitation number. The resulting equation possesses a steady state which solely depends on the dimension of the excitation number subspaces. A comparison to the results obtained from the numerically exact propagation of the Schrödinger equation shows good agreement.

Formulation in the physical degrees of freedom
We consider atoms trapped in a ring-shaped one-dimensional lattice with N sites. We will see, that despite its simplicity this model offers interesting insights into the dynamics of laser-driven Rydberg gases which are usually carried out in a gas cloud. The Hamiltonian is given by where σ (k) x and σ (k) z are the usual Pauli spin matrices acting on the internal state of the atom located at site k and σ (N +1) x . The first part (H 0 ) describes the interaction of the atoms with the laser within the rotating-wave approximation. H int accounts for interaction between Rydberg atoms located different sites. The atomic states are identified as |g k ≡ |↓ k and |r k ≡ |↑ k and 1 + σ (k) z /2 is the projector on the Rydberg state at site k. We are interested here in atomic Rydberg states that interact via the van-der-Waals interaction. In this case we have γ = 6 and the parameter β = C 6 /a 6 determines the interaction strength with C 6 being the van-der-Waals coefficient [27] and a the lattice spacing. The transition |g k → |r k is resonantly (zero detuning) driven with the Rabi frequency Ω. In case of a single atom per site Ω is equal to the single atom Rabi frequency Ω 0 . If each site is occupied by the same number (N 0 ≫ 1) of atoms Ω corresponds to the collective Rabi frequency √ N 0 Ω 0 . In this case the laser couples the product state |g k ≡ [|g k ] 1 ⊗ ... ⊗ [|g k ] N 0 to the superatom ‡ In practice this means that the interaction energy has to be much larger than Ω.
is the symmetrization operator). The derivation and the validity of Hamiltonian (1) is thoroughly discussed in Refs. [10,16].
We are now interested in the regime of the so-called perfect blockade which was introduced in Refs. [13,16]. In this limit we replace the van-der-Waals interaction by an interaction potential whose value is β ≫ Ω for nearest neighbors and 0 for atoms which are further apart. This means that we effectively consider the case in which γ in eq. (1) is infinite. The perfect blockade model is then valid provided that 1 ≫ Ω/β = a 6 Ω/C 6 ≫ 1/64, i.e. Ω has to be much smaller than the nearest neighbor interaction but at the same time much larger than the next nearest neighbor interaction. In practice this means β ≃ 10Ω. For such a choice the further results, which will be obtained under the assumption of the perfect blockade, are virtually indistinguishable from those obtained for van-der-Waals interacting atoms. Such good agreement was also reported in Ref. [13] and is a speciality of the one-dimensional nature of our system; for comparison: In two dimensions the perfect blockade model was valid provided that 1 ≫ Ω/β ≫ 1/8 which is actually difficult to satisfy. The perfect blockade approach has the advantage that certain aspects of the system, e.g. its steady state and derived quantities, can be calculated analytically. Within the perfect blockade the only states |φ which participate in the dynamics are those which satisfy The restriction to this set of states is the manifestation of the strong interaction among the Rydberg atoms. The states |φ can contain any number of Rydberg excitations between 0 and n max , which is the next integer smaller than N/2. The laser Hamiltonian H 0 can only couple states whose excitation number differs by one. Moreover, due to the differences in the spatial distribution of the excitations only certain states are connected by the laser. A convenient way to illustrate the coupling is a graph as shown in fig. 2. Here vertices in the same column contain the same number of Rydberg atoms and adjacent columns are connected by the laser Hamiltonian H 0 . A similar approach to the representation of an interacting many particle system was employed by Altshuler et al. in Ref. [28]. Here a system of interacting fermions was represented as a graph whose vertices are represented by Fock states. Coupling between these states was caused by the interaction. In our case the interaction determines the vertices and the single particle Hamiltonian (H 0 ) determines their coupling. Fig. 2 shows the fully reduced graph for N = 10 where only states which are invariant under reversal (k → N − k + 1) and cyclic shifts (k → k + 1) of the sites of the ring lattice are considered (see Ref. [16]). The state in which all atoms are in the ground state, i.e. |0 = N k=1 |g k , is contained in this fully symmetric set. Solving the Schrödinger equation of Hamiltonian (1) with this (vacuum) state as initial state and calculating the mean number of Rydberg atoms with |Ψ(t) = exp(−iHt) |0 eventually gives rise to the data presented in fig. 1.

Hamiltonian in excitation number space
Due to the strong interaction a formulation of the problem in terms of the physical degrees of freedom (localized atoms) seems disadvantageous. Instead, it appears more natural to base the discussion on the graph depicted in fig. 2. Initially the system is localized on the leftmost vertex of the graph, i.e. it resides in the vacuum state |0 = N k=1 |g k . Once the laser is turned on, coupling to neighboring vertices is established and the propagation of population through the graph sets in. Eventually, the mean Rydberg number is determined by the probability ρ m (t) of the system to be in the subspace containing m excitations: ρ m (t) is hereby the probability density derived from |Ψ(t) and integrated over each column of the graph. Note, that expressions (3) and (4) are equivalent.
In order to make the excitation number spaces explicitly appear in the equations, we use the following matrix representation of the wave function: The In this representation the Hamiltonian becomes where the operators C n−1,n and C † n−1,n connect the subspaces which contain n − 1 and n Rydberg atoms. Each of the subspaces contains dim n states. Hence, C n,n+1 is a dim n × dim n+1 -matrix. The block structure of the Hamiltonian (7) is a direct consequence of the property of the laser Hamiltonian to couple only subspaces whose excitation number differs by one.
The projection operator onto the space containing n Rydberg atoms is given by |n n|, and we can define Ψ n = n | Ψ . This allows us to rewrite Hamiltonian (7) more compactly as |m m| H |n n| = nmax−1 n=0 C n,n+1 |n n + 1| + C † n,n+1 |n + 1 n| .
3. Time-evolution in the excitation number subspace

Time-evolution of the projection operators
Our goal is to calculate the time-evolution of the projection operators |m m|. This will eventually enable us to calculate the quantities ρ m (t). To this end we consider the Heisenberg equation of motion which can be formally integrated to yield For small times τ we obtain up to second order With |m m| 0 ≡ |m m| the first commutator evaluates to where we have used C m−1,m = C † m,m−1 . The second double-commutator becomes The first commutator contains terms of the form Here we have exploited that C † m−1,m is a real matrix. The interaction between the Rydberg atoms manifests itself in the structure of the matrices C † m−1,m which were initially constructed in a product basis in which the single atoms constitute the fundamental degrees of freedom. The strong interaction, however, favors collective excitations which are complex superpositions of the single atom excitations. It is thus reasonable to assume that the single atom degrees of freedom are so strongly mixed that the entries of C † m−1,m can be regarded as uncorrelated. In this case the expression where we have assumed that the complex numbers e i(φ βm −φα m−1 ) are randomly (uniformly) distributed. This allows us to employ the relation | N k=1 e iα k | ≈ √ N since in case of randomly distributed α k we are just dealing with a random walk in two dimensions. The constant c m,m−1 relates to the mean value of the entries of the matrix C m,m−1 .
The same line of argument holds true for terms stemming from the double commutator which are of the form Since the entries of the matrices C m,m+1 and C m+1,m+2 are not correlated, their product is again a matrix with randomly distributed elements. The magnitude of these terms can be estimated by employing again the picture of the random walk in two dimensions. The modulus of (19) then approximately evaluates to c m,m+1,m+2 √ dim m dim m+2 /dim, where c n,m+1,m+2 is a constant related to the mean value of the entries of C m,m+1 C m+1,m+2 . Qualitatively different, however, are the terms of the form where the matrix C m,m+1 appears twice. Since the matrix elements of C m,m+1 are uncorrelated only the diagonal elements of the matrix product yield on average a nonzero value, hence Moreover, since the results cannot depend on the choice of the basis functions spanning a given m-excitation subspace, we can say that [κ m,m+1 ] αm = κ m,m+1 and hence Estimating κ m,m+1 ∝ dim m+1 the modulus of these (diagonal) terms is proportional to dim m dim m+1 /dim. We now return to eq. (15). We neglect all terms of the form (16) and (19) keep only the dominant ones, i.e. those which contain products of the form C m,m+1 C † m,m+1 . By this we obtain Here we have used from which follows that We can thus make the ansatz κ m,m+1 = u m dim m+1 and κ m+1,m = u m dim m and find, omitting the t = 0 argument of ρ m (0), This equation neglects coherent processes which have effectively been eliminated by the neglect of terms of the form (16) and (19). Thus eq. (25) cannot be valid for arbitrary small values of τ as here certainly coherent effects dominate the evolution. Instead, the interval τ has to be chosen sufficiently large such that terms of the form (21) dominate all other contributions whose importance diminishes due to the summation of complex numbers with random phases. Eq. (25) is thus a map that propagates the vector ρ(t 0 ) by a 'coarse-grained' timestep τ , i.e. ρ(t 0 ) → ρ(t 0 +τ ). τ is thereby chosen much smaller than the typical timescale which governs the evolution of ρ m (t).

The steady state
The steady state is defined through i.e. it is a fix point of the mapping (25). The mapping contains the unknown coefficients u m which contain information about how adjacent excitation subspaces are connected. Fortunately, in order to determine the steady state their knowledge is not necessary. It is only required that u m = 0, which is always the case. The solution of eq. (26) is given by Thus it is only the number of states contained in a given excitation number subspace that determines the steady state. It is actually possible to calculate the dimension of the subspaces with fixed number of excitations, analytically. This is done by counting all states of the single atom product basis that contain m excited atoms and obey condition (2). The result is Note, that this result encompasses all possible basis states and not only the fully symmetric ones which constitute the graph ( fig. (2)). It is interesting to see how ρ steady m behaves in the limit of large N, where we have n max = N/2. It is convenient to introduce the variable α = m/N which is the number of Rydberg atoms divided by the number of sites. Using Stirling's formula we can approximate eq. (28) by This function has a very pronounced peak and, for large N, we can approximate ρ steady m (α) by a Gaussian. The position of the maximum of the function is the solution of the equation where ln(x) is the natural logarithm. Since the number of sites N is taken to be very large, we can neglect the last two terms, and then obtain that the function (29) assumes its maximum at Around this peak value the squared width of ρ steady m (α) is given by Hence, for large N the probability distribution in particle number space is strongly peaked with an overwhelming weight on α max . The mean number of Rydberg atoms in the steady state is thus expected to ben r = α max × N and the fluctuation should vanish. The value ofn r is slightly bigger than the values reported in Refs. [13,16]. As we will see in the next section this is due to the particular choice of the initial state.
The distribution ρ steady m (α) contains the full statistics of the Rydberg atom number count. Strong interactions are known to have an effect on the counting statistics leading to a sub-Poissonian distribution [29,30] of the Rydberg number. A measure for this is given by the Mandel Q-parameter which is negative/positive for a sub-/super-Poissonian distribution of the Rydberg atom number count. In the steady state we find which shows the expected sub-Poissonian behavior. For the sake of completeness let us consider the case of non-interacting atoms. Here one obtains for the probability density and hence the probability density in excitation number space performs an oscillatory motion at all times.

Numerical results
We are now going to compare the results that we obtained from the previous section to the actual data obtained from a numerical propagation of the Schrödinger equation. In order to make the numerical solution feasible we massively exploit the symmetry properties of the system. In all numerical calculations we refer to a set of basis states which are invariant under cyclic shifts and reversal of the lattice sites as has been used in Ref. [16]. I.e. we operate only in a subspace of the space which is spanned by all states being compatible with the perfect blockade condition (2).

Evolution into the steady state
Let us start by inspecting the evolution of ρ n when choosing the vacuum as initial state, i.e. ρ 0 = 1, and N = 25. The result is shown in fig. 3a. For t ≤ 5 we observe a welldefined wave packet which propagates through the excitation number space performing an oscillatory motion. For longer times the amplitude of the oscillations decreases, however, the wave packet remains localized. There are still significant fluctuations visible. In particular in the interval 15 ≤ t ≤ 25 remnants of the initial oscillations can be observed. In comparison to the non-interacting case which is shown in fig. 3b the localization in excitation number space is apparent. As expected from eq. (35) here the wave packet exhibits coherent oscillations with large amplitude.

The steady state and its dependence on the initial condition
We now proceed by monitoring ρ n (t) over a time-interval in which the number of Rydberg atoms shows the steady state behavior -here we choose 100 ≤ t ≤ 104. The  result is depicted in fig. 4 for two values of N, 20 and 25. The thin green curves show individual snapshots of ρ n (t) taken at different times. In addition we present also the average of ρ n (t) taken over the considered time-interval. The fluctuations around this average decrease significantly with increasing N. This is in accordance with the behavior of the Rydberg number n r whose fluctuations around the mean value also diminish as N increases (see Ref. [16]). This supports the assumptions that in the limit of very large N indeed a steady state with extremely little fluctuations is established. For both values of N shown in fig. 4 a comparison to the steady state result (27) (thick red curve) reveals a shift of the probability distribution to smaller n. These deviations appear to stem from the particular choice of the initial state: The state |Ψ(0) = |0 is localized at the leftmost vertex of the network. In this region of the graph the matrices C n,n+1 are, however, not actually random since the 'randomness' is caused by the interaction which has little or no effect when the number of Rydberg atoms is only very small, i.e. n = 0, 1, 2. That this 'edge effect' appears to be indeed the cause of the deviation of the probability distribution from ρ steady n is corroborated by the data shown in fig. 5. Here we present the same plot as in fig. 4 but the initial state has to be chosen from the subspace containing 5 excitations, e.g. it is located in the central region of the graph. The effect is not only a much better agreement of the data with ρ steady n but also a significant decrease of the fluctuations about the average of ρ n (t). This behavior is generic for initial states chosen from excitation number subspaces with large dimensions.
Large deviations of the steady state from eq. (27) are again encountered when the initial state is located close to the right hand edge of the graph, i.e. when its number of Rydberg atoms is close to n max (see fig. 5). Hence, it becomes evident that eq. (25) is not unconditionally valid. It constitutes a reliable approximation only if the initial state belongs to a particle number subspace of sufficiently large dimension. That it also works well, when the vacuum is chosen as the initial state is not evident.

Reduced density matrix in excitation number space
So far we have only studied the probability density distribution in excitation number space. Further insights can be gained by examining the reduced density matrix in excitation number space ρ ex as this quantity eventually determines the outcome of a measurement of the number of Rydberg atoms. Its elements are defined by In fig. 6 we present a snapshot of |ρ ex | for a system of 25 sites at the time t = 140. The initial state was |0 . We clearly observe that the entries of the main diagonal dominate the off-diagonal entries. Hence, there is negligible coherence between excitation number subspaces which have a large dimension. Consequently, the reduced density matrix is well approximated by a classical mixture ρ ex ≈ nmax n=0 ρ n |n n|. The strong interaction between the atoms in conjunction with the laser driving erases the phase relation between excitation number subspaces. So tracing out all degrees of freedom but those being relevant for the measurement of the Rydberg number, leaves us with a density matrix of a completely mixed state. This gives actually the impression that the steady state we observe is a state with maximal entropy, since only the dimension of the excitation number subspaces determines the outcome of the measurement (see eq. (27)). In fact we are dealing with a pure state at all times and only the particular measurement we are performing gives us the impression of observing a completely mixed state.

Connection with the microcanonical ensemble
Note that all the information about the steady state could as well have been obtained by considering a microcanonical ensemble. Here, the microstates are just given by the zero energy eigenstates {φ} of the interaction Hamiltonian H int defined through condition (2)) and the steady state given in eq. (28) can be obtained directly by counting the number of these microstates. The fundamental assumption underlying the microcanonical ensemble is that each of the microstates has equal weight. This assumption can clearly not be justified in the absence of the laser, even though all states have the same energy. Only when the laser is present, the microstates {φ} defined by (2)) become strongly mixed and are thus no longer eigenstates of the system. However, they no longer possess strictly zero energy but are rather distributed over an energy interval which is centered at zero and whose width is proportional to NΩ. In other words, although the laser produces a widening of the energy window occupied by the states, it also provides the necessary ingredient that eventually allows the system to thermalize, the equiprobability of the microstates.
At this point we want to remark that the microcanonical prediction to the steady state (28) involves all zero energy eigenstates of H int . On the other hand, the numerical results of section 5 involve only a subset of all accessible states, i.e. the fully symmetric set of eigenstates. The observed agreement between the two approaches suggests that the number of fully symmetric eigenstates with a given excitation number m is proportional to dim m .

Summary and conclusions
We have investigated the origin of the steady state value of the Rydberg number which is exhibited in a laser driven Rydberg gas after an initial transient period. Starting from Heisenberg's equation we have derived an effective equation of motion for the probability density in excitation number space. This effective equation of motion which is coarse-grained in time exhibits a steady state. When comparing this steady state to actual numerical simulations excellent agreement is found provided that the initial state was chosen from an excitation subspace with sufficiently large dimension. In case of an initial state containing a very small/large number of excitations still a steady state is established, however, deviations from the analytical result are obtained.
We have visualized the system by a graph whose vertices are represented by eigenstates of the interatomic interaction. Coupling between the vertices is established by the laser-atom interaction. A similar mapping was applied in Refs. [28,31] where interacting fermions were studied by means of a graph. Here, a transition between localized and delocalized eigenstates has been found to take place as a function of the interaction strength. In our system we are in the regime of strong interaction and the eigenstates are delocalized throughout the entire graph. The observed localization in excitation number space and hence also the observation of a steady state value of n r is a purely statistical effect, owed to the strongly peaked function dim m .
It would be interesting to see whether a distribution of the Rydberg number count similar to eq. (27) and -more specifically -the calculated values for the mean Rydberg number and the Mandel Q-parameter can be observed in actual experiments. Studying the shape and the temporal evolution of the distribution should yield insights into how this steady state is established as the interaction strength increases. Since our simple model is not expected to be valid in higher dimension experiments with Rydberg atoms in lattices could also help here to clarify whether and, if so, how a steady state is established.