Persistent non-Gaussian correlations in out-of-equilibrium Rydberg atom arrays

,


I. INTRODUCTION
Free field theories describe the dynamics of fields with no interactions.These are known as Gaussian as their path integral description contains some form of Gaussian function [1,2].In the presence of interactions, however, the system's ground state often develops strong non-Gaussian correlations.While Gaussian states are highly structured and can be understood using a variety of theoretical techniques, it is the non-Gaussian states that often play a key role as resources for universal quantum computation and enhancing the efficiency of a range of quantum information protocols, including quantum teleporation, communication, sensing, metrology and quantum error correction [3][4][5][6][7][8][9][10][11][12][13][14][15].
Recent works have investigated the behavior of quantum states under non-equilibrium dynamics, sparked by the intriguing question: What happens to an interacting state when interactions suddenly vanish, e.g., after a global quench of the system?It has been shown that generic closed systems, governed by quadratic Hamiltonians, swiftly relax to Gaussian states, regardless of their initial condition in contrast to recent results established in finely tuned open quantum system [16,17].This can be viewed as an example of a "quantum central limit" theorem [18,19].To explain the emergence of Gaussianity, several mechanisms have been proposed, such as spatial scrambling and canonical transmutation, the latter suggesting that Gaussian components of the initial system act as a Gaussian bath, suppressing non-Gaussianity [16,20].These mechanisms have been used to describe the decay of non-Gaussianity in recent experiments on 87 Rb superfluids trapped in a double well potential [21,22].While these studies have provided crucial insights into the process of relaxation in quantum manybody systems [23], they are restricted to systems with effectively non-interacting degrees of freedom, which do not exhibit "full" thermalization but only relax towards a Generalized Gibbs Ensemble [24].It is thus important to understand the role of Gaussianity in interacting systems, which can exhibit chaotic dynamics and thermalization.In particular, it is important to understand if and how non-Gaussianity could be protected in such many-body systems when they are taken out of their equilibrium state.
In this paper, we show that Rydberg atom arrays [25] provide a versatile experimental platform for realizing and manipulating non-Gaussian correlations far from equilibrium.We show that quenching the atoms between different ordered phases allows to explore two very different regimes of correlations.On the one hand, our setup allows to observe how non-Gaussian correlations build up We focus on the Z2 and Z3 ordered phases and quenches between them: blue arrow U Z 2→3 indicates a quench from Z2 ordered phase into the Z3 phase, while the red arrow represents the reverse quench U Z 3→2 .(c) The two types of quenches lead to strikingly different dynamical behaviors.During the U Z 2→3 quench, the system is initially free and Gaussianity grows until it reaches typical values in a random state.By contrast, during the U Z 3→2 quench, the state remains strongly interacting and it is pinned to a highly non-Gaussian manifold.We will quantify this picture in Secs.III and IV using precise Gaussianity measures.
as the system undergoes thermalizing dynamics from an initial, nearly free-fermion state.On the other hand, it is possible to "lock" the system in a strongly non-Gaussian state, which evades both Gaussification and thermalization at late times.This intriguing non-Gaussian regime is found to be remarkably robust, e.g., even to quenching the system across a quantum phase transition.Our proposal can be readily implemented in Rydberg atom experiments [26][27][28][29], which have recently realized the required types of ordered states and protocols for probing correlations in out-of-equilibrium dynamics.
The remainder of this paper is organized as follows.In Sec.II we introduce our model of Rydberg atom arrays and summarize our main results, which are illustrated in Fig. 1.In Sec.III we introduce our approach for quantifying Gaussianity, using two complementary diagnostics: the Wick decomposition in terms of local correlations that are experimentally accessible, as well as a more general measure based on reduced density matrix and variational optimization [30].We apply these tools to obtain the Gaussianity phase diagram of a one-dimensional Rydberg atom array.In Sec.IV we study dynamics of Gaussianity under global quenches connecting different regions of the phase diagram, whose underlying mechanism is eluciated in Sec.V.In Sec.VI we demonstrate the resilience of our results in the face of potential experimental errors, such as local impurity potentials and longer-range interactions.Our conclusions are presented in Sec.VII, while Appendixes contain technical details on the choice of operators in Wick's decomposition, finitesize scaling analysis and role of the boundary conditions.

II. THE MODEL
We consider a one-dimensional (1D) periodic chain containing N Rydberg atoms.Each atom is modelled as a two-level system, where |0⟩ represents an atom in the ground state and |1⟩ is an excited (Rydberg) state.The atomic array is governed by the Hamiltonian [31]: where is the local density operator, and P=|0⟩⟨0| is the projector on the ground state at a given site.The flipping between the ground and excited states is described by the Rabi frequency Ω, U is the chemical potential (detuning), and V is the next-nearest neighbor interaction between atoms excited to Rydberg states, see Fig. 1(a).Unless specified otherwise, due to the periodic boundary conditions, we will restrict our calculations to the zero momentum sector.The model in Eq. ( 1) is applicable to the strong Rydberg blockade regime, where the magnitude of the nearest-neighbor van der Waals interaction is much larger than all other couplings of the model.[25].The strong Rydberg blockade imposes the kinetic constraint n i n i+1 = 0, which forbids Rydberg excitations on adjacent sites.This constraint is enforced in the Hamiltonian by dressing the σ x i with projectors P = |0⟩⟨0| on the neighboring atoms.This prevents the Rabi flip term from generating nearest-neighbor excitations, such that states • • • 11 • • • with neighboring atoms simultaneously excited are projected out of the Hilbert space.The constraint results in the effective model in Eq. ( 1), which is also known as the PXP model [32,33].
The interplay of U and V terms gives rise to a rich phase diagram sketched in Fig. 1(b).The phase diagram was mapped out with high precision in the numerical simulations in Refs.[34][35][36] and explored in experiments [27,28].Large negative values of the chemical potential U favor excitations on every other site (due to the Rydberg blockade, this is largest density of excitations allowed).On the other hand, positive V value assigns a repulsive potential on the next-nearest neighbors and it favors excitations on every third site.Conversely, large negative V favors excitations on every other site.Thus, the model in Eq. ( 1) hosts two ordered phases represented by states in which Rydberg excitations occupy every second or third site, respectively.Both of these phases are destroyed by sufficiently large positive U , which drives the system into a disordered phase, see Fig. 1(b) [31,37,38].
For the subsequent calculations, unless specified otherwise, we set Ω = 1 and concentrate on the Gaussianity and entanglement properties of the initial state as we quench the Hamiltonian between Z 2 , Z 3 ordered phases, Fig. 1(b).We will show that the choice of the initial state and realization of the quench can have dramatically different influence on the Gaussianity, as illustrated in Fig. 1(c).For the quench initialized in the Z 2 phase, indicated by U Z2→3 in Fig. 1(b), the Gaussianity, as precisely defined in Sec.III below, is initially low because the pre-quench state can be approximately expressed as a free-fermion state.After the quench, the state becomes progressively more correlated, with its Gaussianity approaching that of a random vector in the same Hilbert space at late times.This behavior is consistent with thermalization dynamics.In contrast, the ground state in the Z 3 phase cannot be expressed as a free-fermion state and hence it has high non-Gaussianity.Moreover, following the quench, the state remains strongly interacting, which occurs due to a lack of thermalization in this case.It is important to note that in both cases, the quench Hamiltonian, regardless of the nature of the ground state, is an interacting, non-integrable Hamiltonian -further contrasting the two regimes.In the following sections, we introduce several metrics of non-Gaussianity and quantitatively support the phase diagram and dynamical be-havior sketched in Fig. 1.

III. GAUSSIANITY PHASE DIAGRAM
A conventional approach for quantifying the Gaussianity of quantum states relies on Wick's theorem [1].This theorem allows to reduce the evaluation of n-point correlation functions in terms of "contractions" (i.e., vacuum expectation values) of pairs of creation and annihilation operators.For any free-fermion system, the Wick's identity for four-point correlators takes the form where the expectation value is with respect to the ground state ("vacuum").One possible definition of Gaussianity W is the extent to which Eq. ( 4) is violated, i.e., the absolute value of the difference between its left-hand and right-hand side.For Gaussian states, we have W = 0.The operators Â, . .., D are understood to be singlesite fermionic creation and annihilation operators fi , f † j , which obey the anticommutation relation { fi , f † j } = δ ij .As our model in Eq. ( 1) is expressed in terms of spin variables, it will be convenient to work with spin operators rather than fermionic ones, which can be accomplished by applying the Jordan-Wigner transformation.
In order to distinguish the Gaussianity between Z 2 and Z 3 ordered phases, we choose Â = f † 1 , B = f1 , Ĉ = f † 2 , D = f3 , resulting in the following measure of the Wick's decomposition violation: where σ ± j ≡ (σ x j ∓ iσ y j )/2 are the standard spin raising and lowering operators at site j, and ρ denotes the ground state of the system (which could be either a pure state or a density matrix).This particular choice of operators Â, . .., D is justified in Appendix A; it will reveal the difference between the Z 2 phase, where W ≈ 0 across the entire phase, and the Z 3 phase where the deviation from Wick's decomposition is of order unity, W ∼ O(1).
The ambiguity in the choice of operators in the Wick decomposition can be eliminated by employing variational optimization techniques to measure the minimum distance between the reduced density matrix of a given state and the set of all density matrices associated with free-fermion models [30,[39][40][41].This quantity, dubbed the "interaction distance" D F [30], allows for a more general characterization of the state's Gaussianity by quantifying its deviation from the closest free-fermion model defined in an arbitrary basis.Interaction distance is a property of the reduced density matrix ρ describing the subsystem A of a bipartite system A ∪ B. For the total system in a pure state |ψ⟩, the reduced density matrix ρ = tr B |ψ⟩⟨ψ| is obtained by tracing out the subsystem B, and its eigenvalues ρ k define the so-called "entanglement spectrum", E k = − ln ρ k [42].Using the entanglement spectrum, the interaction distance is defined as where l ϵ l is the entanglement spectrum of a free-fermion system, given in terms of single-fermion modes ϵ l and their occupations n [43].The sum runs over the many-body entanglement spectrum.The minimization is over the single-particle energies {ϵ}, whose number typically scales linearly with the number of atoms N .
It is worth noting that the entanglement spectrum is naturally dependent on the choice of bipartition.In this work, we consider a bipartition of the system into two equal parts.Our results, however, are not sensitive to the particular choice of partition, as long as both subsystems are of comparable sizes.If one subsystem is much smaller than the other, the number of Schmidt coefficients will be significantly reduced and we expect non-universal behavior of D F .Interestingly, the interaction distance can also be probed with respect to the eigenspectrum of the system -thus probing the thermal properties of a given Hamiltonian.This analysis can be done and reveals that in both ordered regimes, the Hamiltonian is interacting.Generally, however, this is more cumbersome to perform with increasing system size due to the exponential scaling and does not reveal the distinct differences in the ground state between the two regimes.
Intuitively, D F represents the minimum distance between the reduced density matrix of a given quantum state ρ and the density matrix of the closest free-fermion model defined on the subsystem A. Crucially, the freefermion model is defined up to an arbitrary unitary transformation on A, which makes D F basis independent.This allows to quantify the Gaussianity of a quantum state without the need to search for suitable operators in W as done in Eq. ( 5).However, W has the advantage over D F in that it is expressed in terms of local correlations that are amenable to experimental measurements.Thus, Eq. ( 5) provides a more practical way of detecting non-Gaussianity in the lab.
The Gaussianity phase diagram for the ground state of the Hamiltonian in Eq. ( 1) is presented in Fig. 2 for a range of U and V values.The phase diagram was obtained using both the interaction distance D F and the Wick's theorem violation W, shown in panels (a)-(b).In panel (c) we also show the von Neumann entanglement entropy, S(ρ) = − k ρ k ln ρ k .All the quantities were computed for the reduced density matrix corresponding to the subsystem A being one half of the chain.Fig. 2 reveals excellent qualitative agreement between all three metrics, in particular between interaction distance and Wick's decomposition.The phase boundaries are in good agreement with Refs.[34,36], suggesting weak finite-size effects.For large and negative chemical potential U , there are two competing ordered phases, Z 2 and Z 3 .In particular, for large and positive values of V , the ground state is the Z 3 ordered state.The quantum phase transition from Z 3 to Z 2 ordered state occurs at around |V | ∼ −U/3.In between these two ordered phases, we expect a narrow intermediate commensurate phase [28,31,34], which is difficult to resolve in small systems used in Fig. 2, but this phase will be unimportant for our discussion.
Figures 2(a)-(b) reveal a stark contrast between the two ordered phases in terms of the Gaussian nature of their ground-state correlations: while the Z 2 ground state is approximately a non-interacting Gaussian state with both D F and W close to zero, the Z 3 ground state is nearly maximally interacting, non-Gaussian state.The notion of "maximally interacting" can be made precise by noting that D F , as a trace distance between density matrices, has an upper bound, which has been conjectured to be D max [44].In the Z 3 phase in Fig. 2, D F attains values very close to this upper bound, strongly suggesting that it is not possible to express the Z 3 ground state as a Gaussian state of free-fermionic modes.Finally, we note that the entanglement entropy in Fig. 2(c) also captures some features of the phase diagram, but it does not sharply distinguish between the Z 2 and Z 3 phases.Thus, the interaction distance and local Wick decomposition are essential to gain a complete understanding of non-Gaussianity, both in equilibrium as well as out-of-equilibrium, as we show next.

IV. PERSISTENT NON-GAUSSIAN CORRELATIONS UNDER QUENCH
Previously, we have seen that the two competing ordered phases, Z 2 and Z 3 , are the extreme points on the Gaussianity spectrum: while the Z 2 ground state represents a nearly-free fermion state, the Z 3 state is maximally interacting.It is natural to inquire about the temporal evolution of Gaussianity following a sudden  quench between these phases.According to the standard scenario of thermalization in a closed system [45], under quench dynamics, particularly across a quantum phase transition, the system should lose memory of its initial state and equilibrate towards a maximally entangled state.To test this expectation, we study the spreading of entanglement and non-Gaussian correlations when the system is prepared in the ground state of the Hamiltonian (1) for some value V ≡ V i .We then quench the Hamiltonian to some different value of V ≡ V f ̸ = V i .By varying V i and V f we can access different ordered states and post-quench Hamiltonians.For simplicity, we keep U the same in the initial and post-quench Hamiltonian and postpone the discussion of its role to Sec.VI.
Figure 3(a) contrasts the growth of entanglement entropy for the U Z2→3 quench vs. U Z3→2 quench.In the first case, the system exhibits thermalization, as confirmed by the fast growth of entropy towards its saturation value when it reaches the thermal state.A key indication of thermalization is the volume-law scaling behavior of the saturation value of entanglement entropy, S ∞ ∝ N , consistent with Fig. 3(a).In contrast, quenches from the Z 3 state lead to non-thermalizing dynamics, as seen in the strongly suppressed growth of entropy in Fig. 3(a).
In Fig. 3(b) we illustrate how Gaussianity changes in time when we prepare the system in an approximately Gaussian Z 2 state at t = 0.In particular, when the postquench Hamiltonian is in the Z 3 phase (e.g., V f = 8), the deviation from Gaussianity sharply increases from zero and quickly reaches the saturation value of D ∞ F ≈ 0.03, which coincides with interaction distance of a random vector [44].This is consistent with thermalizing dynamics at infinite temperature in Fig. 3(a), where the state at late times becomes similar to a random vector.Note that this scenario is very different from Ref. [22], where the initial state was chosen to be non-Gaussian, but the Hamiltonian itself is quadratic and induces the development of Gaussian correlations over time.
Conversely, for the non-Gaussian Z 3 initial state in Fig. 3(c), we see that the previous scenario does not hold.In this case, there is persistent non-Gaussianity after the quench, with no sign of decay of the correlations due to interactions.Consequently, the time-evolved state remains highly interacting over the course of quantum dynamics.We note that these results hold for larger system sizes via finite-size scaling and with open boundary conditions -see Appendix B.
We now analyze the U Z3→2 quench in more depth.In terms of spectral properties, we find that there is a single energy eigenstate |E 1 ⟩ of the Hamiltonian that has very high overlap with the Z 3 state, | ⟨E 1 |Z 3 ⟩ | 2 ≈ 1, see Fig. 4(a).The energy of this eigenstate also exactly matches that of the ground state energy of the initial Hamiltonian pre-quench.In fact, the energy of the state is independent of deforming V in the quench Hamiltonian as shown in Fig. 4(b).This implies that the Z 3 state is effectively close to being an eigenstate of the Hamiltonian.This behavior is somewhat reminiscent of quantum many-body scarring [46][47][48], with the exception that here the high overlap exists only for a single eigenstate.Furthermore, numerically we find that the special eigenstate |E 1 ⟩ has lower entanglement entropy than the majority of the spectrum, whereas its interaction distance attains a nearly maximum value.By also plotting the overlap between the Z3 state and the eigenstates of Hamiltonians with V = 8 and V = 1, we see the single eigenstate remains dominant at constant energy, regardless of the value of V , simply transitioning from being an initial ground state to a mid-spectrum state.
Previously, in Fig. 3(c), we saw that non-Gaussianity remains robust for the quench U Z3→2 .In order to experimentally access this behavior, one can study tempo- Overlap between the initial Z3 state (as defined in Eq.( 2)) and energy eigenstates |E⟩ of the Hamiltonian in Eq. ( 1), plotted as a function of energy.Data is for system size N =18 and U = − 15, for three V values given in the legend.In all the cases, the Z3 state has high support (overlap ≈ 1) on a single eigenstate at roughly the same energy.(b) Energy expectation value, E = ⟨ψ|H|ψ⟩, for product states |ψ⟩ = |Z3⟩ and |ψ⟩ = |Z2⟩, plotted as a function of V with fixed U = − 15.We see that |Z3⟩ remains at constant energy for any V , while |Z2⟩ scales linearly.(c) The power spectrum of the correlation function ⟨σ z i σ z i+1 ⟩ evaluated in the time evolved state (raw data shown in the inset).The vertical lines in the plot represent the energy gaps ωij between the states with largest overlap in (a).The gaps align precisely with the peaks in the power spectrum.
ral behavior of local correlation functions, as frequently done in modern ultracold atom experiments [49].For example, the correlation function ⟨σ z i σ z i+1 ⟩, computed in Fig 4(c), reveals persistent oscillations.The characteristic frequencies of these oscillations correspond to the energy differences of the eigenstates with dominant overlap with the initial state of the system.This can be characterized more precisely by the power spectrum [50,51], computed in Fig. 4(c).The dominant frequencies are given by ω 1j = |E 1 − E j | (and their differences), where |E j ⟩, j = 2, 3, . .., denote eigenstates with subleading overlaps with the Z 3 state.Similar oscillations and frequencies can be observed in the quantity W defined in Eq. ( 5), other two-point local correlations, and even in the entanglement entropy.
A simple heuristic argument that gives an approximate value of the oscillation frequency in the limit U, V ≫ Ω can be stated as follows.For the Hamiltonian with U =−15 and V =−5, the ground state is approximately Z 2 product state with energy E Z2 GS ≈ (U + V )N/2.The energy of the Z 3 state is E Z3 ≈ U N/3.For general values of U and V there are no other states with the same energy as Z 3 .For special ratios of U/V , a resonance may occur and other states could have the same energy as Z 3 ; we can prevent this by assuming U and V to be irrational numbers.The oscillations seen in Fig. 4(c) are between |Z 3 ⟩ and states where one of the excitations is moved by a single unit, i.e., the states |101000100100 . ..⟩, |100101000100 . ..⟩ etc., which all have a single 101 pattern.The energy of these states is U N/3 + V , so they are lower in energy by approximately −|V | compared to Z 3 .This predicts that the oscillation frequency is set by V , i.e., the energy difference between Z 3 states and these states containing |. . .101 . ..⟩ in the chain.Thus the energy differences between second and first energy levels are determined by |V |, i.e., ω 12 = |E 1 − E 2 | ≈ |V |, as can be seen in the power spectrum in Fig. 4(c).

V. ORIGIN OF PERSISTENT NON-GAUSSIAN CORRELATIONS
The origin of robust non-Gaussianity associated with the Z 3 state can be more readily understood by considering the evolution in eigenspace overlap for different V presented in Fig. 4(a).We consider the difference between the initial Hamiltonian, H i , and post-quench Hamiltonian, H f .Recalling Fig. 1(b), we restricted to the case where the quench only changes the value of V .Therefore H f = H i + ∆V H nn where H nn = i n i n i+2 .Now consider a quench from the Z 3 state in its respective regime such that there is only occupancy of every third site, therefore H nn |Z 3 ⟩ = 0, irrespective of ∆V .Thus, Hence, upon deforming ∆V , |Z 3 ⟩ remains an eigenstate of H f with the same energy.This may not necessarily still be the ground state and instead may be shifted up in the energy spectrum for sufficiently large ∆V .This means that upon quenching, the initial state remains the same over long periods of time due to its proximity to an eigenstate.As interaction distance is defined only with respect to a given state ρ, it is clear why it does not significantly change over time, despite quenching the system across criticality.This interpretation is supported by the high overlap with a single eigenstate of the quench Hamiltonian in Fig. 4(a).
A similar argument can be made for why the quench from the initial Z 2 state with H f in the Z 3 phase leads to scrambling and thermalization dynamics.As Z 2 has an occupancy on every 2 lattice sites, it "feels" the deformation of ∆V : Considering the form of H i such that the initial state is |Z 2 ⟩, the terms that result in this being approximately the ground state are These competing negative and positive terms mean that, overall, −(|U |+|V |)+∆V may not be much greater than Ω and |Z 2 ⟩ may no longer be approximately an eigenstate like before.This results in the possible scrambling of the initial state into a non-Gaussian state over time.
The further substantiate the previous argument, the underlying mechanism for persistent non-Gaussian correlations can be inferred by considering an effective quench Hamiltonian with five-body interactions.As we are interested in the quench dynamics going from Z 3 into the Z 2 ordered phase, we define the effective Hamiltonian in a regime where U is negative and large, favoring particle creation on all sites.Similarly, we require V to be large and negative as well.Under these conditions, the quench Hamiltonian is given by with |U |≫1, |V |≫1 and we still have n i n i+1 =0.Following a similar procedure as in Ref. [52], we move the quench Hamiltonian into an interaction picture with respect to the next-nearest neighbor term by applying the transformation W † H q W , where Ignoring the rapidly oscillating phases for |V | ≫ 1, we reach an effective Hamiltonian In the largest fully connected sector, the presence of Rydberg excitations on the nearest and next-nearest neighboring sites is prohibited.The effective Hamiltonian corresponds to the PPXPP model with a chemical potential.The Z 2 state still exists as the overall ground state but instead within a small disconnected sector due to the new blockade condition.Meanwhile, due to a large negative U in the effective Hamiltonian, the ground state of the largest sector where the blockade remains respected is Z 3 .Thus, the quench Hamiltonian does not induce delocalizing dynamics when the system is initialized in the Z 3 state and such states are protected against both Gaussification and thermalization.We can therefore conclude that the persistent non-Gaussianity of the Z 3 initial state equivalently arises from the effective blockade mechanism up to the next-nearest neighbor excitations in the interaction picture.Agreement in dynamics was also tested and found numerically between the exact Hamiltonian and effective Hamiltonian.This further supports the notion that the initial state remains approximately an eigenstate of the quench Hamiltonian.

VI. EXPERIMENTAL IMPLICATIONS
With the possible quantum information applications, it is important to test the robustness of the non-Gaussification against external perturbations.This is particularly important because our results rely on the quantum superpositions of states with degenerate energies in the Z 3 and Z 2 phases.External perturbations may result in the superposition collapsing into an energetically favorable product state, thus removing any non-Gaussian correlations.Here we focus on three types of effects that are relevant for experimental implementations: (i) the stability against a single site magnetic field or impurity εn i with magnitude ε; (ii) the effect of changing the chemical potential U during the quench; (iii) the effect of long-range van der Waals interactions that are present in real systems but were neglected in Eq. (1).Fig. 5 shows the results obtained when we add the impurity term εn i to the Hamiltonian in Eq. ( 1) on site i = 4.We choose this site along the chain as it is found to have the most substantial effect on the results providing a qualitative lower-bound in robustness.Despite the presence of an impurity, we see that qualitative features of the phase diagram remain preserved with an impurity strength ε=10 −4 , see Fig. 5(a).This is a magnitude of error much larger than the detuning resolution of current quantum technology [53].Furthermore, perturbations are generally characterised by their proportionality to the ground state gap.We find this order of magnitude to be comparable to the energy gap of the system (which decreases with N ).This demonstrates that the non-Gaussian characteristics are protected nearly up to the same order as the energy gap in the system.It is natural that any larger magnitude of error would disrupt this as one would no longer be probing ground-state physics.Taking a single point in this diagram, marked by the cross, we find the ground state still possesses high overlap with the superposition state Z 3 .Consequently, the non-Gaussian correlations persist when quenching in the Z 2 phase (with error still present in the quench Hamiltonian), as seen in Fig. 5(b).For this point, the non-Gaussianity remains robust for impurity strengths up to ε∼10 −3 .Furthermore, in Figs.5(c)-(d) we test the robustness of D F and W for this point with varying impurity strength and system size.The non-Gaussianity is seen to be more pronounced in smaller system sizes.
While our presented analysis assumed that only V is changed during the quench, we have numerically verified that the non-Gaussian correlations also remain robust upon simultaneous changes in U .This can be understood via the following argument.Consider modulating both V and U , then . This is instead equivalent to quenching horizontally in the phase diagram in Fig. 2; therefore, if ∆U is such that one remains in the regime where |Z 3 ⟩ is approximately the ground state, the state remains an eigenstate and the non-Gaussianity remains robust.This is illustrated by the black dashed line in Fig. 5 where in changing V during the quench, we also take ∆U = 5.On the other hand, if ∆U is large enough to transition from the Z 3 regime, thermalization occurs.This stability against small changes in U makes the non-Gaussianity effect robust against possible experimental imperfections.
Finally, our idealized model in Eq. ( 1) neglects the long-range van der Waals forces that are invariably present in real systems of Rydberg atoms [27][28][29].Thus, it is important to verify our conclusions still hold in the full model describing the Rydberg atom experiments [35]: Note that, in contrast to Eq. ( 1), here we keep the factor 1/2 in the Rabi term and set V =1 in order to facilitate comparison with the literature.In Fig. 6(a), we first recompute the Gaussianity phase diagram of the long-range model with relevant parameters taken from the experimental papers [28,29].Similar to the truncated model in Eq. ( 1), the full model also realizes both Z 3 and Z 2 phases.The phase diagram in Fig. 6(a) is in good agreement with that given in Ref. [35].We then prepare the state in one phase (indicated by a red cross) and perform a quench into the other phase.As illustrated in Figs.6(b)-(c), the results are consistent with those of the UV model, where the Z 3 state preserves its non-Gaussian correlations.For the Z 2 initial state, the thermalization time scale is longer than in Fig. 3 due to the smallness of the energy gap in the chosen units for the Hamiltonian (10).Taking this a step further, we introduce experimental error into our calculation.Ref. [53] states that there are approximate errors of ≈ 0.1µm in the spatial position of sites along the Rydberg chain.We can factor this into our by modulating i, j → i+δ i , j+δ j in Eq. 10 (so, numerically, δ is randomly sampled from a normal distribution between ±0.02).We find that the results still hold well when taking the initial ground state from a disordered Hamiltonian with only a slight decrease in D F as shown by the orange line in Fig. 6(b).More-so, if one assumes the perfect Z 3 state can still be prepared, we find the perfect results still hold irrelevant of the disordered quench Hamiltonian -adding a degree of robustness as it demonstrates the error only factors into the initial state preparation.Overall, the main features of our results are present in the full Rydberg model, suggesting that persistent non-Gaussianity could be observed with the existing experimental technology [28,29].

VII. CONCLUSIONS
Quantum states can exhibit a Gaussian or non-Gaussian nature, depending on the degree of interaction between the system's constituent parts.In this work, we have investigated the phenomenon of Gaussification, wherein non-Gaussian states undergo transformation into Gaussian states during quench dynamics in quantum many-body systems.We have shown that Ryd-berg atom arrays provide a versatile platform where this behavior can be probed with available experimental techniques.Perhaps more intriguingly, we have demonstrated that the Rydberg blockade, intrinsic to such systems, gives rise to states with remarkably robust non-Gaussian correlations that persist far from equilibrium, e.g., as the system is quenched across the quantum phase transition.We have elucidated the origin of this behavior by analyzing quenches between Z 2 and Z 3 phases, which exhibit either scrambling dynamics or suppression of thermalization due to the effective Rydberg blockade mechanism.We formulated a criterion based on Wick's decomposition that incorporates local correlations, providing a practical method for observing (non-)Gaussianity in experiment.This finding was further corroborated via variational optimization and computing the minimal distance between the reduced density matrix of the ground states belonging to different ordered phases and the set of all free-fermion density matrices defined on the same subsystem.
Our results highlight the richness of quantum state complexity in systems evolving under constrained dynamics, and they provide three contributions to the broader quantum information framework.Firstly, our findings show the existence of robust non-Gaussian states in Rydberg systems, which are well-known resources for enhancing quantum information protocols.Recent studies have shown that non-Gaussian states act as magic states [14], facilitating the construction of a universal gate set.Our findings propose a route towards accessing and utilising these robust states in the commonly explored Rydberg systems, moving closer to the realisation of universal quantum devices, even in the presence of long-range interactions and local impurity potentials.Secondly, we demonstrate how Rydberg systems naturally generate and manipulate Z 3 and Z 2 states, which have a stark contrast in terms of their Gaussianity, due to the interplay between detuning and interactions.Particularly, we find that the Z 3 state is maximally non-Gaussian while simultaneously robust against thermalization, which could serve as a qutrit basis for quantum memories [54,55].Lastly, our results provide a constructive example that diverges from the typical Gaussification scenario of Refs.[18,23], thus illustrating the possibility of richer types of dynamical behavior facilitated by the Rydberg blockade.6) and Eq. ( 5), respectively.System size N = 15 and the Wicks decomposition is computed for operators on sites 7, 8, 9 in the bulk of the chain.

1748958.
three to allow for the Z 3 phase.This ensures that excitations will be localized at the edges and the bulk behavior will accurately reflect the underlying physics.Wick's decomposition, implemented using formula (5), can conveniently probe this.Here, due to OBCs, we evaluate local correlations only within the bulk of the system in order to avoid the boundary effects where the Rydberg blockade is slightly weaker.The Gaussianity phase diagrams, depicted in Fig. 7, demonstrate excellent agreement with the PBC results in terms of both interaction distance and local Wick's decomposition in Eq. (5).
Finally, we investigated the interaction distance in larger system sizes undergoing quench dynamics from the Z 3 phase to Z 2 phase.As can be seen in Fig. 8, the non-Gaussianity becomes more pronounced with increasing system size at moderate times t ≲ 15.As the system size increases, the oscillation amplitudes are observed to decrease.Thus, in large system sizes, we anticipate essentially constant D F value with small fluctuations around it.

FIG. 1 .
FIG. 1.(a) Schematic description of a 1D Rydberg atom model.(b) Sketch of the phase diagram in the U -V plane.We focus on the Z2 and Z3 ordered phases and quenches between them: blue arrow U Z 2→3 indicates a quench from Z2 ordered phase into the Z3 phase, while the red arrow represents the reverse quench U Z 3→2 .(c) The two types of quenches lead to strikingly different dynamical behaviors.During the U Z 2→3 quench, the system is initially free and Gaussianity grows until it reaches typical values in a random state.By contrast, during the U Z 3→2 quench, the state remains strongly interacting and it is pinned to a highly non-Gaussian manifold.We will quantify this picture in Secs.III and IV using precise Gaussianity measures.

5 FIG. 2 .
FIG. 2. Gaussianity phase diagram for the ground state of the model in Eq. (1) as a function of U and V .The color scale represents the value of interaction distance DF in (a), the Wick's decomposition violation W in (b), and the entanglement entropy S in (c).All data was obtained by exact diagonalization for N =18 atoms on a ring with periodic boundary conditions.

FIG. 3 .
FIG. 3. Temporal behavior of entanglement and Gaussianity for quenches U Z 2→3 and U Z 3→2 , previously indicated in Fig. 1(b).The chemical potential is held fixed at U =−15.(a) Growth of entanglement entropy for different system sizes.The top three lines represent U Z 2→3 (specifically, Vi=−5 → V f =8), while the bottom three (overlapping) lines are for the reverse U Z 3→2 quench.For the U Z 2→3 quench, the saturation entropy obeys the volume law scaling with system size, indicating thermalization.By contrast, U Z 3→2 quench leads to strongly non-thermalizing dynamics, as evidenced by a complete lack of entropy growth.(b)-(c) Temporal behavior of Gaussianity measured by interaction distance.In (b), we quench from the Z2 ground state (Vi=−5) to a range of V f values spanning both Z3 and Z2 phases.The top plateau value corresponds to the interaction distance of a random state, D random F ≈ 0.03[44], consistent with thermalization observed for U Z 2→3 in (a).(c) Similar to (b) but for Z3 initial state (Vi=8).The persistent large value of DF is consistent with an absence of entanglement spreading for U Z 3→2 quench in (a).Data in panels (b)-(c) is for system size N =18.

FIG. 4 .
FIG.4.The nature of the non-Gaussian quench U Z 3→2 .(a) Overlap between the initial Z3 state (as defined in Eq.(2)) and energy eigenstates |E⟩ of the Hamiltonian in Eq. (1), plotted as a function of energy.Data is for system size N =18 and U = − 15, for three V values given in the legend.In all the cases, the Z3 state has high support (overlap ≈ 1) on a single eigenstate at roughly the same energy.(b) Energy expectation value, E = ⟨ψ|H|ψ⟩, for product states |ψ⟩ = |Z3⟩ and |ψ⟩ = |Z2⟩, plotted as a function of V with fixed U = − 15.We see that |Z3⟩ remains at constant energy for any V , while |Z2⟩ scales linearly.(c) The power spectrum of the correlation function ⟨σ z i σ z i+1 ⟩ evaluated in the time evolved state (raw data shown in the inset).The vertical lines in the plot represent the energy gaps ωij between the states with largest overlap in (a).The gaps align precisely with the peaks in the power spectrum.

N i=1 −
|U |n i −|V |n i n i+2 with U, V ≫ Ω.When then quenching with H f in the Z 3 regime, suitable tuning of positive ∆V then means that given H f , |Z 2 ⟩ is no longer an eigenstate due to the competing factors of U, V and ∆V .More concretely, the final Hamiltonian when acting on |Z 2 ⟩ has a term proportional to −(|U | + |V |) + ∆V .

FIG. 5 .
FIG. 5. Resilience of non-Gaussianity against a local impurity, εn4, added to Eq. (1).(a) The phase diagram obtained using local Wick's decomposition, W(ρ), with impurity strength ε = 10 −4 at system size N = 12.The cross marks the point (U, V )=(−4, 10.5) studied in (b)-(d).(b) The time evolution of the interaction distance DF , when quenching the initial ground state at (U, V )=(−4, 10.5) in the Z3 phase with the Hamiltonian parameters (U, V )=(−4, −6) in the Z2 phase.Both the initial and quench Hamiltonian contain impurity ε.Data is for system size N = 12 with the impurity potentials shown in the legend.Black dashed line is with no error (ε = 0) but instead taking the initial ground state at (U, V )=(−15, 8) and quenching at (U ′ , V ′ )=(−10, −5), which demonstrates persistent non-Gaussinity even with a change in U .(c)-(d) Wick's decomposition and interaction distance, respectively, of the ground state at (U, V )=(−4, 10.5), as a function of impurity strength for several system sizes.The results were computed with exact diagonalization without resolving translation symmetry.

FIG. 6 .
FIG. 6.(a) Gaussianity phase diagram of the long-range Rydberg model in Eq. (10) for system size N =12 and fixed V =1.(b)-(c) Temporal behavior of DF when quenching from the point indicated by the red cross in the Z3 (Z2) phase to the other phase.The different behavior of non-Gaussianity for the two types of quenches, seen in Fig. 3, is reproduced.Blue lines indicate using the pure Rydberg Hamiltonian Eq. 10.Meanwhile the orange line introduces further experimental error by using a Hamiltonian with spatial disorder such that i, j → i + δi, j + δj where δ is a site dependent and taken randomly in the range [−0.02, 0.02].The results are averaged over 100 realisations.Both the pure and disordered Rydberg Hamiltonian show persistent non-Gaussinity.

FIG. 7 .
FIG. 7. Gaussianity phase diagram for open boundary conditions (OBC).The are analogous to Figs. 2(a-b) with interaction distance and Wicks decomposition computed according to Eq. (6) and Eq.(5), respectively.System size N = 15 and the Wicks decomposition is computed for operators on sites 7, 8, 9 in the bulk of the chain.