Ubiquitous quantum scarring does not prevent ergodicity

In a classically chaotic system that is ergodic, any trajectory will be arbitrarily close to any point of the available phase space after a long time, filling it uniformly. Using Born’s rules to connect quantum states with probabilities, one might then expect that all quantum states in the chaotic regime should be uniformly distributed in phase space. This simplified picture was shaken by the discovery of quantum scarring, where some eigenstates are concentrated along unstable periodic orbits. Despite that, it is widely accepted that most eigenstates of chaotic models are indeed ergodic. Our results show instead that all eigenstates of the chaotic Dicke model are actually scarred. They also show that even the most random states of this interacting atom-photon system never occupy more than half of the available phase space. Quantum ergodicity is achievable only as an ensemble property, after temporal averages are performed. It is generally believed that most eigenstates of quantum chaotic models are ergodic. In this work, the authors disprove this by showing that all eigenstates of the Dicke model in the chaotic regime are scarred, and that ergodicity is an ensemble property, achievable only in the temporal average.

A striking feature of the quantum-classical correspondence not recognized in the early days of the quantum theory is the repercussion that measure-zero structures of the classical phase space may have in the quantum domain. A recent example is the effect of unstable fixed points, that cause the exponentially fast scrambling of quantum information in both integrable and chaotic quantum systems [1][2][3][4][5] . Another better known example is the phenomenon of quantum scarring [6][7][8] . As a parameter of a classical system is varied and it transits from a regular to a chaotic regime, periodic orbits that may be present in the phase space change from stable to unstable. These classical unstable periodic orbits can get imprinted in the quantum states as regions of concentrated large amplitudes known as quantum scars. Even though the phase space may be densely filled with unstable periodic orbits, they are still of measure zero, which explains why it took until the works by Gutzwiller 9 for their importance in the quantum chaotic dynamics to be finally recognized.
Quantum scarring was initially observed in the Bunimovich stadium billiard 10 and soon in various other one-body systems [11][12][13] giving rise to a new line of research in the field of quantum chaos 8,[14][15][16][17][18][19][20][21][22] . The recent experimental observation of long-lived oscillations in chains of Rydberg atoms 23 , associated with what is now called "many-body quantum scars", has caused a new wave of fascination with the phenomenon of quantum scarring [24][25][26][27][28] . While the interest in many-body quantum scars lies in their potential as resources to manipulate and store quantum information, a direct relationship between them and possible structures in the classical phase space has not yet been established.
Halfway between one-body and many-body models, one finds systems such as two-dimensional harmonic oscillators and the Dicke model 29 , where quantum scars have also been observed 30,31 . In the first case, the model is not fully chaotic and scarring can be understood as an extension of the regular orbits 32,33 . The Dicke model, on the other hand, has a region of strong chaos, where the Lyapunov exponents are positive and the level statistics agrees with the predictions from random matrix theory 34 . The model describes a large number of two-level atoms that interact collectively with a quantized radiation field and was first introduced to explain the phenomenon of superradiance 35,36 . It has been studied experimentally with cavity assisted Raman transitions 37,38 , trapped ions 39,40 , and superconducting circuits 41 .
In this work, we investigate the intricate relationship between quantum scarring and phase-space localization in the superradiant phase of the Dicke model. Even though both phenomena are often treated on an equal footing, the connection is rather subtle. Scarring refers to structures that resemble periodic orbits in the phase-space distribution of quantum eigenstates, while phase-space localization implies that a state exhibits a low degree of spreading in the phase space. Here, we demonstrate that scarring does not necessarily imply significant phase-space localization.
In systems studied before, scarred eigenstates were thought to be a fraction of the total number of eigenstates. In contrast to that, we show that deep in the chaotic regime of the Dicke model, all eigenstates are scarred. Their phase-space probability distributions always display structures that can be traced back to periodic orbits in the classical limit. Yet, we find eigenstates that are highly localized in phase space and eigenstates that are nearly as much spread out as random states, although none of them, including the random states, can cover more than approximately half of the available phase space.
In addition to the analysis of quantum scarring and phasespace localization, we also provide a definition of quantum ergodicity. This is done using a measure that we introduce to quantify the level of localization of quantum states in the phase space. We say that a quantum state is ergodic if its infinite-time average leads to full delocalization. Under this definition, stationary quantum states are never ergodic, while random states are, and coherent states lie somewhere in between.

Results
Dicke model and chaos. The Hamiltonian of the Dicke model is written asĤ where ℏ = 1. It describes N two-level atoms with atomic transition frequency ω 0 interacting with a single mode of the electromagnetic field with radiation frequency ω. In the equation above, a (â y ) is the bosonic annihilation (creation) operator of the field mode,Ĵ x;y;z ¼ 1 x;y;z are the collective pseudo-spin operators, withσ x;y;z being the Pauli matrices, and γ is the atom-field coupling strength.
The eigenvalues j(j + 1) of the squared total spin operator b J 2 1 J 2 x þĴ 2 y þĴ 2 z specify the different invariant subspaces of the model. We use the symmetric atomic subspace defined by the maximum pseudo-spin j ¼ N =2, which includes the ground state. When the Dicke model reaches the critical value γ c ¼ ffiffiffiffiffiffiffiffi ωω 0 p =2, it goes from a normal phase (γ < γ c ) to a superradiant phase (γ > γ c ). Our studies are done in the superradiant phase, γ = 2γ c , and we choose ω = ω 0 = 1. The rescaled energies are denoted by ϵ = E/j. For the selected parameters, ϵ GS = −2.125 is the ground-state energy.
The classical Hamiltonian, h cl (x) in the coordinates x = (q, p; Q, P), is obtained by calculating the expectation value of the quantum Hamiltonian under the product of bosonic Glauber and pseudo-spin Bloch coherent states x j i ¼ q; p j i Q; P j i (see Methods) and dividing it by j. The effective Planck constant ℏ eff = 1/j 42 determines the resolution of the quantum states on the four dimensional phase space M. We are able to work with large system sizes (j~100 and Hilbert space dimensions D~6 × 10 4 ), due to the use of an efficient basis that guarantees the convergence of a broad range of eigenvalues and eigenstates 43,44 (see Methods). The Dicke model displays regular and chaotic behavior. For the Hamiltonian parameters selected in this work, the system is in the strong-coupling hard-chaos regime for ϵ > −0.8 (see Supplementary Note 1).
Quantum scarring. The (unnormalized) Husimi function of a stateρ is defined as QρðxÞ ¼ x h jρ x j i, which is the expectation value of the density matrix over the coherent state x j i centered at x. This function is used to visualize how the stateρ is distributed in the phase space. Quantum scars are localized around the classical periodic orbits in an energy shell of the phase space. To visualize the scars, we consider the Husimi projection over the classical energy shell at ϵ, By integrating out the bosonic variables (q, p), the remaining function can be compared with the projection of the classical periodic orbits on the plane of atomic variables (Q, P).
Identifying all periodic orbits that generate the scars of a quantum system is extremely challenging. We were able to identify two families of periodic orbits for the Dicke model, which we denote by family A and family B 45 . By calculating the overlap of the eigenstates with tubular phase-space distributions located around these orbits 8 , we selected twelve eigenstatesρ k ¼ E k j i E k h j scarred by those two different families. In Fig. 1 we plot their Husimi projections e Q k ¼ e Q ϵ k ;ρ k at ϵ k = E k /j along with the corresponding periodic orbit of each family. Family A (solid blue line in Fig. 1a1-a6) contains the periodic orbits of lowest period of the Dicke model, which emanate from one of the two normal modes around a stable stationary point at the ground-state energy. Family B (solid red line in Fig. 1b1-b6) arises from the other normal mode around the same point. Scarring is clearly visible in all panels of Fig. 1. The quantum states are highly concentrated around the classical periodic orbits. This happens even in the chaotic region of high excitation energy, where the classical dynamics is ergodic, as seen in Fig. 1a5, a6, b5, b6. Notice that the eigenstates may be scarred by more than one periodic orbit. In fact, as we showed quantitatively in ref. 45 after introducing a measure of scarring, an eigenstate may even be scarred by periodic orbits of different families. This means that different eigenstates exhibit different degrees of scarring.
It is evident from Fig. 1 that the degree of delocalization of the eigenstates in phase space also varies. The Husimi distribution of the eigenstates in Fig. 1a5, a6, for instance, is not entirely confined to the two periodic orbits drawn in blue. This contrasts with the high density concentration that the eigenstate in Fig. 1a4 shows around the plotted unstable periodic orbits. To quantify these differences, we introduce a measure of the degree of localization of a quantum state in the classical energy shell.
Scarring vs. phase-space localization. To measure the localization of a state in a Hilbert-space basis indexed by some letter n, the most commonly used quantity is the participation ratio P R 46-48 . Its inverse is given by P À1 R ¼ P n P 2 n , where P n is the probability of finding the state in the n'th basis vector. By analogy, we introduce a measure of localization in phase space that employs as basis the overcomplete set of coherent states within a single energy shell M ϵ ¼ fx ¼ ðq; p; Q; PÞjh cl ðxÞ ¼ ϵg, so that we replace the sum ∑ k with a three-dimensional surface integral R M ϵ ds over M ϵ . For a given x 2 M ϵ , the probability P x of finding the stateρ in the coherent state x j i is given by the Husimi function QρðxÞ. With these replacements, we finally obtain a measure of phase-space localization Lðϵ;ρÞ given by ds QρðxÞÞ 2 =VðϵÞ is a normalization constant and ds is the volume of M ϵ (see Methods). The measure Lðϵ;ρÞ is an energy-restricted second moment of the Husimi function 49 . It is related to the second-order Rényi-Wehrl entropy 47 , which, in turn, was shown for the Dicke model 50 to be linearly related to the first-order Rényi-Wehrl entropy 51 .
The value of Lðϵ;ρÞ indicates the fraction of the classical energy shell at ϵ that is covered by the stateρ. It varies from its minimum value Lðϵ;ρÞ $ ð2π_ eff Þ 2 =VðϵÞ, which indicates maximum localization, to Lðϵ;ρÞ ¼ 1, which implies complete delocalization over the energy shell. The former occurs for coherent states, and the latter happens if QρðxÞ is a constant for all x 2 M ϵ , in which case the projection e Q ϵ;ρ ðQ; PÞ is also constant for the allowed values of Q and P.
All eigenstates in Fig. 1 have values of L k ¼ Lðϵ k ;ρ k Þ below 1/2. For the eigenstates in Fig. 1a1, a2, a4, b3, these values are very small, since the eigenstates are almost entirely localized around the plotted periodic orbits. The value of L k is larger in Fig. 1a3, because at the center of the diagram, there is an unstable stationary point 2 , which produces a one-point scar in addition to the scar associated to the blue orbit. The localization measure is larger in Fig. 1b1, b2 simply because in these cases the phase space is very small. It is also larger for the states in the high energy region in Fig. 1a5, a6, b4, b5, b6, because they spread beyond the marked periodic orbits. As the energy increases and one approaches the chaotic region, more unstable periodic orbits emerge in the classical limit, enhancing the likelihood that a single quantum eigenstate gets scarred by different periodic orbits. We stress, however, that even for those high-energy states with larger values of L k , the drawn periodic orbits cast a bright green shadow that is clearly visible in the Husimi projections.
It is important to make it clear that scarring and localization, despite related, are not synonyms. Of course, there is no eigenstate with a large value of L k that would at the same time have a high value of our scarring measure, but there is more to the relationship between these two concepts. A highly localized eigenstate is scarred by few periodic orbits of a particular family and therefore has a high value of the scarring measure for that family, although it has low values for the other families of periodic orbits 45 . Furthermore, even the most delocalized eigenstates are still scarred, but now by periodic orbits from different families 22 .
In Fig. 2, we take a step further in the analysis of localization and scarring. In the large panel in Fig. 2a, we show L k against energy for all eigenstates between ϵ GS and ϵ = 0.06. This plot is equivalent to a Peres lattice 52 for expectation values of observables, as used in studies of chaos and thermalization. In the low-energy regular regime, L k is organized along lines that can be classified with quasi-integrals of motion linked with classical periodic orbits 53 . Conversely, as the system enters the chaotic region at higher energies, the distribution of L k becomes dense and looses any order. Notice, however, that all eigenstates in the chaotic region have values of L k much lower than 1, mostly clustering below 1/2.
The value L $ 1=2 marks a limit on the spreading of any pure state in the high energy shells of the phase space. To show this, we build random states R ϵ j i ¼ P k c k E k j i, where c k are complex random numbers from a Gaussian energy profile centered at energy ϵ (see Methods). The values of the localization measure Lðϵ;ρ ϵ Þ for four different random statesρ ϵ ¼ R ϵ j i R ϵ h j centered at increasing energies ϵ between −0.6 and −0.1 are given in Fig. 2r1-r4, and indeed L $ 1=2. This upper bound on the phasespace delocalization is not due to quantum scarring, but to quantum interference effects.
The panels in Fig. 2r1-r4 display the projected Husimi distributions e Q ϵ;ρ ϵ ðQ; PÞ of the four random states and those in Fig. 2s1-s22 Fig. 1 (see also ref. 45 ), we have not identified the periodic orbits associated with Fig. 2s1-s22, but the visible circular patterns are clear evidence of periodic orbits. Their existence is supported by the shape of the Husimi projections and by knowing the generic direction of the classical Hamiltonian flow. The patterns display all the features of periodic orbits: they always cross the line P = 0 perpendicularly, they are symmetric along both the horizontal Q and P axes, and they visibly form closed loops. There is no quantum effect other than scarring that would produce such patterns. One therefore deduces that deep in the chaotic region of the Dicke model, all quantum eigenstates are scarred, although they have different degrees of scarring, as seen by the patterns, and they have different degrees of localization, ranging from strong localization (L $ 0:1) to states that are nearly as much delocalized as random states.
Our results are sharpened as one approaches the semiclassical limit. The patterns indicating scarring do not fade away as the system size increases. Quite the opposite, as j = 1/ℏ eff increases, the periodic orbits get better defined in the Husimi projections (cf. the figures for j = 30 and j = 100 in the Supplementary Note 3). To study the dependence of the level of phase-space localization on system size, we show the distributions of L for random states (Fig. 2b) and for eigenstates in the chaotic region (Fig. 2c) for different values of j. For the random states, L is concentrated around 1/2 and the width of the distribution decreases as j increases, corroborating that L ¼ 1=2 is indeed the delocalization upper bound. In contrast, for the eigenstates in the chaotic regime, the distributions are skewed and broader. The tail at small values of L k does not change as j increases, showing that the highly localized states persist, while the portion of the states with large L k decreases, suggesting that for large system sizes, none of the eigenstates reach values of L k > 1=2.
The ubiquitous scarring revealed by our studies motivates the question of whether scarring in other quantum models is also the rule. We have found hints in the literature suggesting that our findings may actually be quite general. For example, in ref. 22 , the authors reconstruct significant portions of the spectrum of a quantum chaotic system using only periodic orbits. This means that all of these eigenstates are described by those periodic orbits and are therefore scarred. In ref. 20 , the authors claim that the great majority of the eigenstates of the hydrogen atom in a magnetic field may be related to periodic orbits, indicating that scars must be the rule. But to provide a definite answer, a phasespace analysis similar to the one presented here is needed.
Quantum ergodicity. We have so far discussed two conceptsquantum scarring and phase-space localizationthat are related, but are not equal. How about their relationship with quantum ergodicity? In the classical limit, a system is ergodic if the trajectories cover the energy shell homogeneously. We then adopt the same definition for quantum ergodicity. To quantify how much of the energy shell is visited on average by the evolved statê ρðtÞ ¼ e ÀiĤ D tρ e iĤ D t , we consider the infinite-time average 54,55 , and compute Lðϵ;ρÞ Lðϵ; ρÞ with Eq. (3). If the whole energy shell at ϵ is homogeneously visited byρ, then Lðϵ;ρÞ ¼ 1. We thus say that a quantum stateρ is ergodic over the energy shell ϵ if Lðϵ;ρÞ ¼ 1. According to this definition, all stationary states in the chaotic region of the Dicke model are non-ergodic, since Lðϵ k ;ρ k Þ ¼ Lðϵ k ;ρ k Þ ≲ 1=2, as shown above. How about non-stationary states, such as coherent states or random states, are they ergodic? We study the evolution of initial coherent states Ψð0Þ j i¼ x 0 j i ¼ P k c k E k j i with mean energies ϵ = −0.5, that are in the chaotic region (see Methods). We select both coherent states that are highly localized and delocalized in the energy eigenbasis, with the degree of delocalization measured by the participation ratio P R ¼ ð P k jc k j 4 Þ À1 . Their energy distributions are shown in Fig. 3a1-g1. The components of the states with low P R are bunched around specific energy levels (Fig. 3a1, b1), exhibiting the comb-like structure typical of scarred states 6 . As P R increases, the coherent states become more spread in the energy eigenbasis, looking more similar to the random state Ψð0Þ j i¼ R ϵ j i shown in Fig. 3h1, whose mean energy is also ϵ = − 0.5.
The evolution of the survival probability, S P ðtÞ ¼ j Ψð0Þ h je ÀiĤ D t Ψð0Þ j ij 2 , for the coherent states with low P R leads to large revivals before the saturation of the dynamics 56 , as seen in Fig. 3a2, b2, c2. This contrasts with the evolution of the coherent states with large P R , such as those in Fig. 3d2-g2, and the evolution of the random state in Fig. 3h2. In these cases, the approach to the asymptotic value of S P (t) is much smoother and exhibits the so-called correlation hole, which corresponds to the ramp towards saturation. The correlation hole reflects the presence of correlated eigenvalues and is a quantum signature of classical chaos [57][58][59] .
The values of Lðϵ;ρÞ for the states in Fig. 3a1-h1 are indicated in Fig. 3a3-h3. The random state is indeed ergodic, reaching L $ 1. For the coherent states, L increases as P R does, but even for the states with the largest values of the participation ratio, L is still slightly under 1. The analysis of the dependence of L on system size is done in Fig. 3i, j. The distribution of L for random states (Fig. 3j) gets narrower and better centered at L ¼ 1 as j increases, confirming that these states indeed behave ergodically. The distribution for the coherent states (Fig. 3h) is much broader. The center moves towards larger values as j increases, but it is not clear whether there will ever be a significant portion of the initial coherent states with L $ 1.
In Fig. 3a3-h3, we plot the projected time-averaged Husimi distributions e Q ϵ;ρ ðQ; PÞ for the states in Fig. 3a1-h1. Remarkably, even for the coherent states with high P R , which do not exhibit any comb-like structure in their energy distributions and do not show revivals in the evolution of their survival probability, we still see an enhancement around unstable periodic orbits, as revealed by a careful inspection of Fig. 3d3-g3 and in contrast with the absence of any pattern for the random state in Fig. 3h3. This unexpected manifestation of dynamical scarring can only be observed in phase space, having no identifiable signature in the energy distribution of the initial states or in the evolution of the survival probability. Thus, revivals in the long-time quantum dynamics are signs of a scarred initial state, but lack of revivals do not exclude the presence of scarring, it just indicates that if it exists it is at a low degree.

Discussion
The three main concepts investigated and compared in this work were quantum scarring, phase-space localization, and quantum ergodicity. We showed that for the Dicke model, all eigenstates in the chaotic region are scarred, although with different degrees of scarring and different levels of phase-space localization. Evidently, an eigenstate that is strongly scarred is also highly localized in phase space, but a single eigenstate may be scarred by different periodic orbits and reach levels of delocalization almost as high as a random state. We also showed that any pure stateeven without any trace of scarringis localized in phase space, and that ergodicity is an ensemble property, achievable only through temporal averages. Thus, scarring, localization and lack of ergodicity are not synonyms, although connections exist.
The ubiquitous scarring of the eigenstates does not immediately translate into the breaking of quantum ergodicity. All eigenstates are certainly non-ergodic, since they never reach complete delocalization in phase space, but if a non-stationary state visits on average the available phase-space homogeneously, then it is ergodic. Random states, for example, are ergodic. The analysis of initial coherent states showed that some are heavily scarred, resulting in the strong breaking of ergodicity that translates into the usual revivals of the survival probability. More interesting is the subtle behavior of the majority of the initial coherent states, which do not display revivals in the quantum dynamics or the comb-like structure in their energy distributions, but yet show some degree of scarring.
Analyses that focus on the Hilbert space, such as the energy distribution of the initial states or the fluctuations of eigenstate expectation values in Peres lattices and comparisons with thermodynamic averages, that are often done in studies of the eigenstate thermalization hypothesis (ETH), may miss the ubiquitous scarring observed in this work. For this feature to be revealed, one needs to look at the structures of the states in phase space.
Our results for the Dicke model are, of course, important, due to the widespread theoretical interest in this model and the fact that it is employed to describe experiments with trapped ions and superconducting circuits. But the repercussion of the findings discussed here goes beyond the limit of spin-boson systems by raising the question of whether scarring in other quantum systems is also the rule and not the exception. Our work provides the appropriate tools to address this question. The phase-space method that we developed is applicable to any quantum system that has a tractable phase-space. Whether these studies could eventually be extended also to interacting many-body quantum systems, such as interacting spin-1/2 models, will depend on the viability of their semiclassical analysis, and some recent works give reasons for optimism [60][61][62] .

Methods
where N is the eigenvalue of the modified bosonic number operatorÂ yÂ and m 0 ¼ m x is the eigenvalue of the original collective pseudo-spin operatorĴ x . The modified bosonic vacuum states in The Hilbert space dimension of this basis is given by D ¼ ð2j þ 1ÞðN max þ 1Þ, where N max designates a cutoff of the modified bosonic subspace. This basis allows to work with larger values of j by reducing the value of N max required for convergence of the high-energy eigenstates. With j = 100 and N max ¼ 300 (D = 60,501), we are able to get 30,825 converged eigenstates covering the whole energy spectrum up to ϵ = 0.853. Having converged eigenstates in such a high-energy regime would be infeasible with the usual Fock basis for j = 100 44,67-69 .
Husimi projection and localization measure. To compute the Husimi projection given in Eq. (2) and the localization measure given by Eq.
where x = (q, p; Q, P) and f(x) is a non-negative function in the phase space. For the localization measure, note that Lðϵ;ρÞ À1 ¼ 1 By using properties of the δ function, those integrals are reduced to e f ðQ; PÞ ¼ Z p þ p À dp P q ± f ðq ± ; p; Q; PÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Δðϵ; p; Q; PÞ p ; ð8Þ where q ± are the two solutions in q of the second-degree equation h cl (q, p; Q, P) = ϵ, Δðϵ; p; Q; PÞ ¼ ∂h cl ∂q ðq ± ; p; Q; PÞ and p ± are the two solutions in p of the second-degree equation Δ(ϵ, p, Q, P) = 0. Because of the form of the weight 1= ffiffiffi ffi Δ p , the integral given by Eq. (8) may be computed efficiently using a Chebyshev-Gauss quadrature method.
It is worth noting that a quantum state with a Wigner distribution that is constant in an energy shell will lead to a Husimi function that needs not to be constant within the same energy shell. This is because the Husimi distribution is the convolution of the Wigner distribution with the Gaussian Wigner distributions of the coherent states, which have different energy widths due to the geometry of the energy shells in the phase space. This rather marginal effect may be seen in Fig. 3h3, where there is a barely visible weak concentration towards the center of the plot causing L to be slightly under 1. We stress that this effect is not related to quantum scarring. It is just a manifestation of the phase-space geometry in the Husimi distributions.
Random states. The state R ϵ j i ¼ P k c k E k j i is built by sampling random numbers r k > 0 from an exponential distribution λe −λx and random phases θ k from a uniform distribution in [0, 2π). We use where ν(E) is the density of states, ρ(E) is a Gaussian profile of width jσ centered at energy jϵ, and M ensures normalization. This way, R ϵ j ihas a defined energy center ϵ, where the Husimi projection and the localization measure are calculated 59,70 .

Data availability
All the data that support the plots within this paper and other findings of this study are available from the corresponding authors upon request.

Code availability
All the computational codes that were used to generate the data presented in this paper are available from the corresponding authors upon request.