Quantum scarring in a spin-boson system: fundamental families of periodic orbits

As the name indicates, a periodic orbit is a solution for a dynamical system that repeats itself in time. In the regular regime, periodic orbits are stable, while in the chaotic regime, they become unstable. The presence of unstable periodic orbits is directly associated with the phenomenon of quantum scarring, which restricts the degree of delocalization of the eigenstates and leads to revivals in the dynamics. Here, we study the Dicke model in the superradiant phase and identify two sets of fundamental periodic orbits. This experimental atom-photon model is regular at low energies and chaotic at high energies. We study the effects of the periodic orbits in the structure of the eigenstates in both regular and chaotic regimes and obtain their quantized energies. We also introduce a measure to quantify how much scarred an eigenstate gets by each family of periodic orbits and compare the dynamics of initial coherent states close and away from those orbits.


I. INTRODUCTION
The eigenstates of quantum systems that are fully chaotic in the classical limit may not be completely delocalized in the energy shell. This realization [1] came as a surprise, since in the classical domain the trajectories do fill the phase space densely. The presence of unstable periodic orbits (UPOs) plays a central role in the phase-space representation of the eigenstates [2]. Despite having measure zero, UPOs give rise to quantum scars, which are characterized by enhanced probabilities along the phase-space region occupied by those orbits and thus prevent the uniform distribution of the eigenstates [3]. Quantum scars were first studied in one-body systems [2,[4][5][6][7][8][9][10][11][12] and then in two-dimensional harmonic oscillators [13,14] and the Dicke model [15][16][17]. More recently, they have been connected also with some special non-thermal states of many-body quantum systems [18,19], although in this case the analysis of the classical limit is still missing.
The present work focuses on the effects of quantum scars in the Dicke model. This spin-boson model consists of N atoms collectively coupled with a quantized field. It was introduced to understand the superradiance phenomenon in lightmatter systems [15,[20][21][22][23][24] and it was later used in studies of nonequilibrium dynamics [25][26][27][28][29][30][31][32], including the evolution of out-of-time-ordered correlators [33][34][35], and as a paradigm of the ultra-strong coupling regime in several systems [36][37][38][39]. The model can be studied experimentally with cavity assisted Raman transitions [40,41], trapped ions [42,43], and circuit quantum electrodynamics [44]. It has two degrees of freedom and displays both regular and chaotic behavior [23,[45][46][47][48][49][50] depending on the Hamiltonian parameters and excitation energy. Scarring in the Dicke model was studied for a low number of atoms in [16,51], where an algorithm to identify the classical periodic orbits (POs) was implemented. The technique allowed to match the phase-space probability accumulations of the eigenstates with specific POs. The comparison between quantum and classical phase-space distributions was later extended to a large number of atoms in [17].
Recently, in [52] we showed that even in the chaotic high-energy region, all eigenstates of the Dicke model accumulate around portions of the phase-space energy shells of the corresponding eigenenergies, covering at most half of the available phase space. This was done using a phase-space localization measure similar to the one studied in Ref. [53].
In this work we identify the two fundamental families of POs that emanate from the two normal modes of the groundstate configuration and study extensively their influence in the phase-space localization of the eigenstates. The orbits change from stable in the regular low-energy regime of the model to unstable as one approaches the chaotic high-energy region. We introduce a measure that quantifies the degree of scarring of an eigenstate by a specific PO and use it to select the states concentrated around those two families of POs. We also show that the energies of these states follow the Bohr-Sommerfeld quantization rules for both stable and unstable POs.
The second part of the paper is dedicated to the dynamics of the Dicke model in the chaotic regime. The evolution depends on the proximity of the initial state to the POs. In the classical limit, an initial configuration launched exactly at a UPO will remain over it for an infinitely long time. However, the UPOs form a set of measure zero in the phase space, so their existence does not break the ergodic properties of the whole system and initial configurations picked up at random pass arbitrarily close to all other accessible configurations at the same energy [54]. In the quantum domain, the effects of the UPOs are more dramatic. An initial state defined over a phase space region that includes a short-period UPO should exhibit revivals and, after long times, it should be more likely to be found in the vicinity of that same UPO. Therefore, the state explores an effective reduced volume of the phase space [32] and displays what is known as a dynamical scar [55] in its infinite-time average. In contrast, initial states that are away from short-period UPOs should follow a path to equilibrium according to random matrix theory, maximally covering the phase space at long times. Since we have two families of POs, we can choose initial coherent states that are very close to them and verify these expectations.
The paper is organized as follows. In Sec. II, we present the Dicke Hamiltonian, its properties, and its classical limit.
In Sec. III, we identify the families of classical POs emanating from the ground-state configuration and introduce a measure to determine the degree of scarring of the eigenstates. In Sec. IV, we analyze the effects of scarring in the dynamics of initial coherent states close to those POs. Our conclusions are given in Sec. V.

II. DICKE MODEL
The Dicke Hamiltonian [15] is given bŷ where = 1. The model describes a system of N two-level atoms with atomic transition frequency ω 0 interacting with a single mode of the electromagnetic field with radiation frequency ω. The parameter γ controls the atom-field coupling, a (â † ) is the usual bosonic annihilation (creation) operator of the field mode, andĴ x,y,z = 1 2 N k=1σ k x,y,z are the collective pseudo-spin operators withσ x,y,z being the Pauli matrices.
The eigenvalues j(j + 1) of the squared total-spin operator J 2 =Ĵ 2 x +Ĵ 2 y +Ĵ 2 z determine the different invariant subspaces of the model. We use the symmetric atomic subspace defined by the maximum pseudo-spin value j = N /2, which includes the ground state. The HamiltonianĤ D commutes with the parity operatorΠ = e iπΛ , where the operatorΛ =â †â +Ĵ z + j1 has eigenvalues λ = n + m + j, that correspond to the total number of excitations. The number of photons is given by n and the number of excited atoms by m + j, where m is an eigenvalue of the atomic operatorĴ z .
When γ reaches the critical value γ c = √ ωω 0 /2, the Dicke model presents a second-order quantum phase transition [20][21][22][23]. At this point, it goes from a normal phase (γ < γ c ), where the ground state has no photons and all atoms are in their ground state, to a superradiant phase (γ > γ c ), where the ground state has a finite amount of photons and excited atoms.

A. Classical limit of the Dicke Hamiltonian
The classical Hamiltonian defined over the fourdimensional phase space M of the Dicke model, with coordinates x = (q, p; Q, P ), is constructed using Glauber-Bloch coherent states [16,47,48,50,51,56] |x = |q, p ⊗ |Q, P , which are built as a tensor product of Glauber coherent states for the bosonic sector and Bloch coherent states for the pseudo-spin sector where Z 2 = Q 2 +P 2 , |0 is the photon vacuum, and |j, −j is the state with all atoms in their ground state. The raising (lowering) collective pseudo-spin operator,Ĵ + (Ĵ − ), is defined in the usual wayĴ ± =Ĵ x ± iĴ y . Taking the expectation value of the HamiltonianĤ D under the states |x and dividing it by the pseudo-spin j [32], we obtain the classical Hamiltonian We define the rescaled energy corresponding to h cl as which determines an effective Planck constant eff = 1/j [57]. Depending on the Hamiltonian parameters and excitation energies, different dynamical behaviors of the model are identified, ranging from regularity to chaos. As a case study, we choose ω = ω 0 = 1, coupling parameter in the superradiant phase γ = 2γ c , and pseudo-spin value j = 30 (N = 60). For these parameters, the ground-state energy is given by GS = −2.125. The dynamics is regular up to a value of ≈ −1.7 and chaotic for higher energies [50].

III. QUANTUM SCARRING
Quantum scarring is a phenomenon by which some eigenstates of a quantum system get concentrated around the UPOs that appear in the classical limit of the model. To identify the scarred eigenstates of the Dicke model, we must first identify the POs arising in the classical limit.
A. Families of periodic orbits in the classical limit A PO O with period T is a subset of the phase space, such that The most trivial periodic orbit consists of a single stationary point, where x st (t) = x st (0) for all times t. The classical dynamics given by h cl yield several stationary points that may be located by finding the extrema of the Hamiltonian. For our chosen parameters, there are two stationary points located at the ground state energy GS [56]:  We first consider x GS , which is marked in Fig. 1 with a blue arrow in Figs. 1 (a1) and (a2), and with a red arrow in Figs. 1 (b1) and (b2).
By considering small displacements around the stationary point x GS , we obtain two normal frequencies of the system, Ω A GS and Ω B GS , which are given by where A corresponds to the plus sign and B to the minus sign. For our selection of parameters, Ω A GS = 4.008 and Ω B GS = 0.966. Correspondingly, we have two normal periods T A GS = 2π/Ω A GS = 1.568 and T B GS = 2π/Ω B GS = 6.503. Let us focus first on the normal mode of period T A GS with the trivial PO O A GS = {x GS }. By perturbing this stable stationary point and using a monodromy method [58] to guarantee the convergence, we can find a new PO O A with energy = GS + δ and period T A = T A GS + δT [59]. This procedure can be successively repeated to increasing energies , so that a continuous family of periodic orbits O A is obtained all the way to the chaotic regime, where the orbits become unstable (see App. A for details). This family of POs is denoted by A, and it is shown in Figs. 1 (a1) and (a2), where the color of each PO indicates its energy (lighter colors mean larger energies). We repeat the procedure above for the other normal mode around the ground state with T B GS , yielding the POs of family B, which is plotted in Figs. 1 (b1) and (b2). The energy of each PO in the two families identified increases as the POs grow away from x GS . In Fig. 1 (c), we show the period T A of the POs in A (blue line) and T B for the POs in B (red line) as a function of the energy . For a given energy , the period T B is between three and four times larger than T A , a feature to which we come back to when explaining our next results. In Fig. 1 (d), we display the Lyapunov exponents of those POs as a function of energy. The orbits become unstable when the Lyapunov exponents are different from zero. Notice that for family B, the Lyapunov exponent does not grow monotonically with the energy and there are ranges of high energies where the PO can be stable.
The classical Hamiltonian (5) is invariant under the transformation (q, p; Q, P ) → (−q, p; −Q, P ). This is the classical manifestation of the parity conservation present in the quantum system. The invariance means that if O is a PO, we may find another PO by mirroring where O has the same period, energy, and Lyapunov exponent as O. This transformation yields the mirrored families The two families A and B are the ones that emanate from the stationary point x GS .

B. Scarring of energy eigenstates
The four families A, B, A, and B scar many of the quantum eigenstates, as we show in this subsection.

Measure of scarring degree
The Husimi function of a stateρ can be used to visualize how it is distributed in the phase space. Thus, in order to find the eigenstates scarred by those families of POs, we make use of the (unnormalized) Husimi function of a stateρ, where |x is the coherent state centered at x. In the case of a pure stateρ = |ψ ψ|, the function reduces to To quantify the degree of scarring of a given state by a specific PO, we consider the temporal average of the state's Husimi function along that PO. For a stateρ and the PO O with period T , we define the quantity where x(t) ∈ O and any initial point x(0) ∈ O can be used to perform the average. This measure is similar to the tube phase-space projection introduced in Ref. [10], but now taking into account the temporal behavior of the orbit. From the definition, we have that whereρ is a tubular Gaussian distribution around the PO and x(t) is the coherent state centered at the point x(t). See App. B for an illustration of the Husimi function ofρ O . From Eq. (14), we can see that Qρ O is the overlap of stateρ withρ O . To construct a baseline to compare the value of Qρ O with, we consider a totally delocalized statê comprised of all the coherent states within a single energy is the phase-space volume of the energy shell at . The value of the trace tr(ρ ρ O ) gives the overlap between a totally delocalized state and the orbit. By defining we obtain a direct measure of quantum scarring. A value of P(O,ρ) = 1 indicates that the overlap between stateρ and the PO O is equal to that of a totally delocalized state. Values greater than 1 indicate that the stateρ is scarred by the PO O. A value of P(O,ρ) = 2, for example, says that stateρ is twice as likely to be found near the PO as compared to the delocalized stateρ . Values less than 1 signal that stateρ is less likely to be found near the PO than a fully delocalized state. This avoidance of a specific PO may be regarded as an "anti-scarring" effect. We may quantify how much a stateρ is scarred by the POs of the families described in Sec. III A by defining the numbers for any energy . In words, P A measures how muchρ populates the region of the phase space visited by the POs of families A and A at the energy , while P B does the same for the families B and B.

Overlap of the energy eigenstates with classical periodic orbits
For an eigenstateρ k = |E k E k | of the Dicke Hamiltonian H D with scaled eigenenergy k = E k /j, the numbers measure the scarring produced by the orbits of families (A, A) and (B, B), respectively [60]. In Fig. 2, we show the expectation value of the operator j z =Ĵ z /j for each energy eigenstate, leading to an arrangement known as a Peres lattice [47,61]. This is a convenient way to get information about all the eigenstates in a single picture. In the low-energy regime the points are clearly separated and arranged in a lattice-like fashion, but as the energy increases the structure gets destroyed and the number of points becomes much denser, as typical of chaotic systems. In Fig. 2 (a), we color the points according to the degree of scarring of each eigenstate |E k by the family A, that is, according to the value of P A k . The same is done in Fig. 2 (b) for family B using now P B k . In Fig. 2 (a) and Fig. 2 (b), we also plot with blue and red thick lines the average of the classical variable j z = (Q 2 + P 2 )/2 over each PO from family A (blue) and B (red), where T A,B is the period of O A,B . We see in Fig. 2  Using the Bohr-Sommerfeld quantization rule, we semiclassically quantize the energies of the families A and B. Beginning from the ground-state energy Bottom panels A1-B6: The green shades indicate the projected Husimi distribution QE k (Q, P ) for the corresponding eigenstates marked in panels (a) and (b) (darker colors indicate higher values). The gray curved outlines mark the border of the energy shell. The colored lines in the panels A1-B6 draw the POs: We plot the obtained energies Fig. 2 (a) [Fig. 2 (b)]. These vertical lines coincide perfectly with the individual eigenstates scarred by the POs of the families A and B in the low-energy region, where the POs are stable. As the energy increases and the POs become unstable, we find clusters of scarred eigenstates dis-tributed around the semiclassical energies. The width of these clusters is given by the Lyapunov exponent of the PO at the respective energy, as indicated with horizontal bars shown at the bottom of Fig. 2 (a) and Fig. 2 (b). The width of these horizontal bars is twice the value of the Lyapunov exponents of (red) multiplied by eff = 1/j. This result is in agreement with the Gutzwiller trace formula for the density of states [62] when applied to the POs of the families A and B.
The distribution of scarred states in the spectrum is governed by two numbers: (i) the periods of the POs determine the semiclassically quantized energies E i around which clusters of scarred states appear, and (ii) the Lyapunov exponents control the width of these clusters. These two numbers behave differently for families A and B as explained below.
(i) The periods of the POs in family A are between three and four times shorter than those of B. This has two consequences, P A k is higher than P B k and there is a smaller number of eigenstates scarred by family A as compared to B. The former occurs because the POs from family A are shorter in the phase space than those from family B, thus having smaller overlaps with a delocalized state. Since the denominator in P(O,ρ) in Eq. (17) is the overlap of the PO with a delocalized state, this results in the values of P A k being between three to four times higher than P B k . The latter occurs because a shorter period translates in greater energy separations E i+1 − E i . Thus, there are approximately four times more red vertical lines in Fig. 2 (b), marking E B i , than blue vertical lines in Fig. 2 (a), marking E A i . Because scarred eigenstates appear close to these energies, this results in more states scarred by family B than states scarred by family A.
(ii) The Lyapunov times of the POs in family A are always much larger than their periods, which reflects the fact that the energy differences E A i+1 − E A i are larger than the Lyapunov exponents multiplied by eff . For family B, the period of its POs in the chaotic region are only slightly shorter than the respective Lyapunov times. This allows us to clearly identify the clusters of eigenstates scarred by family A around the semiclassical energies E A i in Fig. 2 (a), but makes it harder to distinguish the ones scarred by B around E B i in Fig. 2 (b), because neighboring clusters may overlap.
We close Sec. III B 2 with the interesting observation that POs from both families may affect the same eigenstate. In Fig. 2 (a) [Fig. 2 (b)], the blue circles outlined by red [red circles outlined by blue] mark the eigenstates where both P A k and P B k are greater than 4. These five states are located at energies where the semiclassical quantizations of both families coincide.

Projected Husimi distribution
We now use the Husimi distribution to visually confirm that the energy eigenstates with large values of P A k and P B k are scarred by the classical POs of families A and B.
For a stateρ and energy , we consider the projection of the Husimi function Qρ over the classical energy shell h cl (x) = by integrating out the bosonic variables (q, p), where x = (q, p; Q, P ). For an eigenstateρ k = |E k E k |, the projection Q k = Q k ,ρ k yields a function depending only on the variables (Q, P ), which can be compared with the projection of the POs over the same plane Q-P (see App. C for details on the computation of this projection).
We select 12 energy eigenstates [52], marked as A1-A6 and B1-B6 in Fig. 2 (a) and Figs. 2 (b), which in addition to having high values of P A,B k lie close to the classical average of j z given by Eq. (20). We plot the Husimi projection Q k (Q, P ) for each one of these eigenstates at the 12 bottom panels of . Scarring is clearly visible in all panels. The quantum states A1-A6 and B1-B6 are highly concentrated around the classical periodic orbits. This happens even in the chaotic region of high excitation energy, as seen for A5, A6, B5, and B6 with k > −0.5, where the classical dynamics is ergodic Interesting features are revealed by the juxtaposition of the Husimi projections and the periodic orbits. For example, the eigenstate B3 shows a significant concentration of probability towards the center at Q = P = 0. This is because close to the energy of this state, specifically at q = p = Q = P = 0 and = −1, there is an unstable stationary point and so the dynamics of the periodic orbit slows down around it. This gets reflected in the eigenstate by the localization of the Husimi distribution in the same region [see Fig. 6 (b) in App. B for an illustration of this effect]. Also noticeable is the fact that the probability distributions of the eigenstates A5 and A6 are not entirely confined to the periodic orbit, but extend beyond it, which contrasts with the high density concentration of the eigenstates A1-A4. This results in lower values of P A k for A5 and A6 as compared to A1-A4.

IV. SCARRING AND DYNAMICS
In this section, we show the effects that the scarred states have over the dynamical properties of non-stationary states. We consider three initial Glauber-Bloch coherent states that have energy in the chaotic regime. One is centered in a point of a UPO of family A, one in a point of a UPO of family B, and the third one is away from the POs of the identified families A, B, A, and B .

A. Initial coherent states
We expand the selected coherent states in the Hamiltonian eigenbasis The amplitudes c k and the degree of scarring of the corresponding eigenstates |E k determine the dynamical properties of the initial states.

(ai)-(aiii) [Figs. (bi)-(biii)] also display vertical lines that mark the semiclassically quantized energies E
given by Eq. (21). We observe the following features for the three initial coherent states: (i) For the initial coherent state |x i located in the PO of family A, the largest components |c k | 2 , indicated as i1-i6 in Fig. 3 (ai), correspond to the eigenstates with large values of P A k and therefore scarred by the POs of family A. Due to the high participation of these eigenstates, we say that the initial coherent state |x i is itself scarred by family A. As visible in Fig. 3 (ai), the LDoS of |x i exhibits a clear comb-like pattern, which is typical of scarred states [3]. This comb-like structure is inherited from the density of states. As seen in Fig. 2 (a), the scarred eigenstates cluster around the semiclassical energies E A i . Because these scarred eigenstates give the largest contributions to the initial state, that is they lead to the biggest components|c k | 2 , the LDoS attains higher values around the energies E A i . One sees that the contributions to |x i from eigenstates non-scarred by the POs of family A are erratically distributed, with small and medium values of |c k | 2 , as evident from the red points in Fig. 3 (bi), which mark the eigenstates according to their values of P B k .
(ii) A similar picture emerges for the initial coherent state |x ii located in the PO of family B. Its largest components |c k | 2 , indicated as ii1-ii6 in Fig. 3 (bii), correspond to the eigenstates with large values of P B k and thus scarred by the POs of family B. The contributions from eigenstates non-scarred by this family are smaller and their values fluctuate randomly. The initial coherent state |x ii is therefore scarred by family B. Its comb-like structure is somewhat visible, but, because the separations E B i − E B i+1 are of the order of the Lyapunov exponents of the POs in family B, it is harder to distinguish it.
As discussed in the previous section, there are more eigenstates scarred by family B than by family A, due to the difference in the periods of the POs between the two families. This difference between the two families is evident if we compare the number of blue circles in Fig. 3 (ai) with the number of red circles in Fig. 3 (bii). The larger number of contributing states from family B explains why the biggest components in Fig. 3 (bii) are smaller than the ones in Fig. 3 (ai).
(iii) For the initial coherent state |x iii , which is located far away from the POs of families A, B, A and B, the largest components |c k | 2 , indicated as iii1-iii6 in Figs. 3 (aiii) and (biii), are not colored by any of the two families, as expected.
In Figs. 3 (ci)-(ciii), we plot the same Peres lattice of the pseudo-spin operatorĵ z that we showed in Fig. 2, but we now use tones of green to indicate the values of |c k | 2 for the coherent states |x i [ Fig. 3 (ci)], |x ii [ Fig. 3 (cii)], and |x iii [ Fig. 3 (ciii)]. The panels make it evident that the components of the three initial states, all covering a wide range of energies within the chaotic regime, are actually localized in different regions of the lattice. The components of |x i (|x ii ) concentrate towards the top (bottom) of the lattice, because this is where the eigenstates scarred by family A (B) tend to be found, while the components of |x iii spread more homogeneously over the middle part of the lattice, as seen in Fig. 3 (ciii). (a) By first intersecting the Husimi distribution of the eigenstate with the energy shell at the respective eigenenergy h cl (x) = k , and then integrating over (q,p), as we did in Fig. 2 A1-B6 [see Eq. (22)] and in [52]. This choice is shown in green in the top of each panel of Fig. 4. (b) By directly integrating over the bosonic variables (q, p) of the Husimi function, dq dpQ k (x), as done in [51,53,64]. This alternative is plotted in orange in the bottom of each panel of Fig. 4.
As expected, the Husimi functions of the eigenstates i1-i6 [ii1-ii6] with the largest components in the coherent state |x i [|x ii ] concentrate along the POs of family A [B] at the corresponding energy. This is more evident with the Husimi functions intersected by the energy shells (top green plots in Fig.  4) than in the complete projected Husimi functions (bottom orange plots in Fig. 4), which are more blurred, because they include all energy shells. This difference shows the advantage of the projection method (a), which makes it easier to distinguish the outline of the POs. Method (a) was first developed in [52] and more details are given in App. C.
Similarly to the eigenstates i1-i6 and ii1-ii6, the Husimi functions of the eigenstates iii1-iii6 that contribute the most to the initial state |x iii are not smoothly distributed either. As seen in Figs. 4 (iii1)-(iii2), the visible concentrations in parts of the phase space suggest that these eigenstates are also scarred by one or more families of POs , although they are different from the families A and B.

B. Survival probability of coherent states
To study the dynamics of the three selected initial states, we consider the survival probability defined as whereÛ (t) = e −iĤ D t is the unitary evolution operator. S P (t) measures the probability of the evolved stateÛ (t) |x to return to its initial state |x for any given time. By expanding the coherent state in the Hamiltonian eigenbasis as in Eq. (23), we can rewrite the survival probability as which can also be written in an integral form, showing that the survival probability is the squared norm of the Fourier transform of the LDoS, as defined in Eq. (25). In general, the evolution of scarred initial states, where the LDoS is fragmented, produces partial revivals at times before the Lyapunov time [1] and the saturation of the dynamics happens at values larger than those for non-scarred states [10,32]. The reason for this behavior is, of course, the low number of eigenstates that contribute to the dynamics.
The survival probability of the initial coherent states |x i , |x i , and |x i are respectively shown in Figs. 5 (ai)-(bi), Figs. 5 (aii)-(bii), and Figs. 5 (aiii)-(biii). The orange curves are numerical results and the red ones are running averages. The initial decay of S P (t) for the three cases is Gaussian, since the envelope of the LDoS is Gaussian, and it then reaches values close to zero, as explained in [32]. The subsequent behavior differs among the three states.
For state |x i , the first revival of the survival probability occurs at t = 2.8 [ Fig. 5 (bi)], which is precisely the average of the periods T A within the energy range of state |x i , T A i . A second revival is observed at t = 2T A i , while the next ones do not follow this period, because T A changes strongly through the energy range of state |x i destroying the periodicity. This may be seen in the left inset of Fig. 1 (c), where the little blue square marks the value of T A i . For longer times, the survival probability equilibrates, fluctuating around its asymptotic value drawn with a horizontal black dotted line [32].
The first revival of S P (t) for state |x ii appears at t = T B ii = 9.5 [ Fig. 5 (bii)], which corresponds to the average period of the POs in the family B within the energy range of state |x ii . The subsequent revivals happens at multiples of this period. The periods T B do not vary as much as T A around = −0.5, as seen in the right inset of Fig. 1 (c), where the value T B ii is marked with a little red square. The slope of T A around T A i is large, while T B actually attains a local minimum close to T B ii . As a result, the revivals of the survival probability of state |x ii follow the period T B ii for much longer than in the case of state |x i .
A great advantage of having identified the families A and B is that we know why and where the revivals should happen. Having access to the number of contributing eigenstates and knowing why there are more contributing states from family B than from family A help us understand also why the survival probability for the initial state |x ii saturates at a lower point than S P (t) for the state |x i .
In contrast to states |x i and |x ii , the initial coherent state |x iii does not show revivals. Instead, the survival probability shows a behavior similar to what we obtain when considering an initial state where the coefficients c k are random numbers. This latter case is indicated with a green line in Fig. 5 (aiii) and Fig. 5 (biii) and it is well described using random matrix theory [29,32]. This is puzzling at first sight, since Figs. 4 (iii1)-(iii6) suggest that the most contributing eigenstates to |x iii , that is eigenstates iii1-iii6 in Figs. 3 (aiii) and (biii), are also scarred. What we have come to understand from the analysis of various initial states is that the onset of revivals depends on two factors: the eigenstates with the largest participations in the LDoS should be scarred by the same family of POs, so that the periods of the POs generating the scars are similar; and these periods should be small in comparison to the Lyapunov times. For a given family of POs, if the Lyapunov time is shorter than the period, the revival does not have time to develop before saturation. This happens if either the PO is very unstable or if it has a very long period. It may therefore be that the eigenstates iii1-iii6 either do not belong to the same family and thus the corresponding POs have very dissimilar periods, or that these periods are so long that the revivals are unable to manifest. These are, however, open questions that can only be answered with the identification of the families associated with those states.

C. Dynamical scarring
We finally analyze how the scarring manifests in the infinite-time average of the coherent states [52] With the functions P A and P B defined by Eq. (18), we may calculate Notice that P A i = 4.77 and P B ii = 2.46 are particularly large, while the rest of the numbers are less or approximately equal to 1. This means that, at any time t, stateÛ (t) |x i [Û (t) |x ii ] is more likely to be found in the vicinity of the orbits of families A and A [B and B] at energy = −0.5 than what one would expect for a state that is completely delocalized in the energy shell. This effect is known as dynamical scarring [55].
We can visually see the dynamical scars by calculating the projected Husimi distributions Q ,ρ [Eq.  Fig. 5 (ciii), which suggest that state |x iii is scarred, but by POs from families other than those identified in this work.
The fact that the survival probability of state |x iii is so well described by random matrix theory [green line in Fig. 5 (aiii) and Fig. 5 (biii)] suggests that even though it has a minor degree of scarring, according to Fig. 5 (ciii), this is not identifiable in the structure of the LDoS. To better understand this feature, we use the ratio between the asymptotic value of the survival probability S ∞ P of |x iii and the asymtpotic value of ). In panels (aiii) and (biii), the green curve represents the analytical expression obtained for the evolution of a random state with a Gaussian energy profile [32]. Panels The ratio R measures the similarity between the LDoS of the coherent state and that of a random state, (see Ref. [32] for an extensive discussion of R). A ratio R = 1 indicates that these distributions are very similar. Lower values imply that there are periodic structures in the LDoS of the chosen coherent states, which are, evidently, absent in the LDoS of a random state [29]. For states |x i and |x ii , R i = 0.41 and R ii = 0.81, respectively, while R iii = 1.01 for state |x iii . These numbers show that the scarring can be traced back to a comb-like structure in the LDoS of states |x i and |x iii , but not of state |x iii . This is consistent with Ref. [52], where we found concentrations resembling dynamical scars in the phase-space projections of even the most random-like coherent states at energy = −0.5, which display no revivals, no comb-like structures in their LDoS, and a path to equilibrium well-described by random matrix theory.
The dynamical behavior of the survival probability of the coherent states can be described by two competing effects. Periodic structures in the LDoS give rise to revivals, while random-like spreading within the Gaussian envelope of the LDoS prevents revivals. In between the two cases, one may find coherent states whose LDoS display a slightly larger participation of some eigenstates separated by a nearly constant energy difference, but this periodic structure does not stand out significantly over the Gaussian envelope, so revivals are not observed. Yet, when one performs infinite-time averages, the scars corresponding to the high participating eigenstates become visible as in Fig. 5 (ciii). A complete description of this kind of coherent states requires the challenging task of identifying the family, or set of families, of POs that scar those high participating eigenstates, so that one can compare their Lyapunov times and periods. This is an interesting idea for a future work.

V. CONCLUSIONS
We identified the two fundamental families of classical periodic orbits (POs) that emanate from the ground state of the Dicke model in the superradiant phase and extensively explored their effects in the quantum domain.
• By introducing a measure of scarring based on the temporal average of the Husimi function of an eigenstate, we were able to identify which eigenstates are scarred by which family.
• By projecting the Husimi functions over the energy shell corresponding to the energy of the eigenstate, we found an effective way to visualize the concentrations around the POs of those families.
• We also showed that knowledge of the periods and Lyapunov exponents of the POs in the two identified families allows us to characterize the distribution of scarred eigenstates. The energies of the eigenstates scarred by the two families cluster around values obtained by means of the Bohr-Sommerfeld quantization rule of the periods, and the width of these clusters is directly related to the Lyapunov exponents of the POs generating the scars.
The dynamical consequences of the presence of eigenstates scarred by the two identified families were studied by considering the survival probability of three representative initial coherent states, two localized in the vicinity of the POs of the two families and one away from them.
• For the two states close to the POs, our detailed knowledge of the families allowed us to understand the revivals, their periods, and the saturation values of the survival probability. By performing the infinite-time averages of the density matrices of these initial coherent states, dynamical scars were observed.
• The third initial state shows a survival probability and a local density of states akin to those of random initial states, even though the dynamics has contributions from eigenstates scarred by families not identified here. The fact that the potential families associated with these scars are unknown to us prevents us from making conclusive statements about this state.
An important extension of the present study would be the systematic identification of more families of POs and the analysis of how they influence the spectrum and dynamics of the model.

ACKNOWLEDGMENTS
We thank D. Wisniacki for his valuable comments and acknowledge the support of the Computation Center -ICN, in particular of Enrique Palacios, Luciano Díaz, and Eduardo Murrieta. SP-C, DV and JGH acknowledge financial support from the DGAPA-UNAM project IN104020, and SL-H from the Mexican CONACyT project CB2015-01/255702. LFS was supported by the NSF grant No. DMR-1936006.
Appendix A: Algorithm to find families of periodic orbits emanating from stationary stable points Given a stable stationary point x GS ∈ M with energy GS and normal period T GS , we iteratively find a continuous family of POs O starting from O = {x GS }. The existence of these families is guaranteed by theorem 2.1 of Ref. [59].
First, we detail a variant of an algorithm known as the monodromy method [58,65]. This is a Newton-Raphson-type algorithm that converges towards a PO given an initial guess for an initial condition and period. This type of algorithms has been extensively studied in several systems. See, for example, Refs. [16,66]. Then, we detail an algorithm to iteratively construct the guesses required to find the POs.
1. The monodromy method: converging to a periodic orbit given a good initial guess Assume we have guessesx andŤ for an initial condition and period of a PO, respectively. We want x =x + ∆x and T =Ť + ∆T , the initial condition and period of an orbit in the same energy shell asx. Denote by Φ x (t) the fundamental matrix associated to the Hamiltonian system h cl and ϕ t (x) the Hamiltonian flow satisfying ϕ t (x) = x(t) (See Also, we can approximate the energy constriction h cl (x) = h cl (x) to first order to get ∇h cl x · ∆x = 0, and, finally, the constriction to stay in the same Poincaré section of constant P as ξ · ∆x = 0, where ξ = (q = 0, p = 0; Q = 0, P = 1) . This last constriction eliminates movement along the flow and increases the stability of the algorithm. The linear constrictions (A2), (A3) and (A4), may be written in matrix form as This overdetermined system of linear equations may be approximately solved by least squares using Moore-Penrose pseudoinversion. The solution for ∆x and ∆T is not exact, but the process may be iterated with new guessesx + ∆x anď T + ∆T that will converge to x and T if the initial guesses were good enough.

Finding families of periodic orbits
We now define an iterative process to find the families of POs. For the first step, we start with a stable stationary point x GS ∈ M with energy GS = h cl (x GS ) and normal period T GS .
Given O , we now explain how to find O with = δ + close to . Pick x ∈ O and define a perturbation δx = a∇h cl (x), where a is a scalar such that h cl (x + δx) = + δ . For the first step, ∇h cl (x GS ) = 0, in which case one may select any direction (we use the q direction), and the stability of the initial stationary point will guarantee that the algorithm converges. Setx where the prev , T prev are the energy and period of the closest previously calculated orbit. This way,Ť is linear extrapolation based on the behavior of the previous orbits.
Using the monodromy method detailed in the previous subsection, we may correct the guessesx andŤ to obtain actual solutions x and T so that the desired PO is Although energy is constrained in the monodromy method we used, this is only to first order, so = h cl (x ) may not be exactly equal to δ + . This is easily fixed by performing additional iterations with smaller |δ | which converge to the desired energy.