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

The classical Hamiltonian defined over the four-dimensional phase space $\mathcal{M}$ of the Dicke model, with coordinates x = (q, p; Q, P), is constructed using Glauber–Bloch coherent states [20, 51, 52, 54, 55, 59]

$$\begin{equation}\left\vert \boldsymbol{x}\right\rangle =\left\vert q,p\right\rangle \otimes \left\vert Q,P\right\rangle ,\end{equation} \tag{ 2 }$$

which are built as a tensor product of Glauber coherent states for the bosonic sector

$$\begin{equation}\vert q,p\rangle ={\mathrm{e}}^{-\left(j/4\right)\left({q}^{2}+{p}^{2}\right)}\enspace {\mathrm{e}}^{\left[\sqrt{j/2}\left(q+ip\right)\right]{\hat{a}}^{{\dagger}}}\vert 0\rangle ,\end{equation} \tag{ 3 }$$

and Bloch coherent states for the pseudo-spin sector

$$\begin{equation}\vert Q,P\rangle ={\left(1-\frac{{Q}^{2}+{P}^{2}}{4}\right)}^{j}\enspace {\mathrm{e}}^{\left[\left(Q+iP\right)/\sqrt{4-\left({Q}^{2}+{P}^{2}\right)}\right]{\hat{J}}_{+}}\vert j,-j\rangle ,\end{equation} \tag{ 4 }$$

where |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, ${\hat{J}}_{+}$ (${\hat{J}}_{-}$), is defined in the usual way ${\hat{J}}_{{\pm}}={\hat{J}}_{x}{\pm}i{\hat{J}}_{y}$.

Taking the expectation value of the Hamiltonian ${\hat{H}}_{\mathrm{D}}$ under the states $\left\vert \boldsymbol{x}\right\rangle$ and dividing it by the pseudo-spin j [36], we obtain the classical Hamiltonian

$$\begin{equation}{h}_{\text{cl}}\left(\boldsymbol{x}\right)\equiv \frac{\langle \boldsymbol{x}\vert {\hat{H}}_{\mathrm{D}}\vert \boldsymbol{x}\rangle }{j}=\frac{\omega }{2}\left({q}^{2}+{p}^{2}\right)+\frac{{\omega }_{0}}{2}\left({Q}^{2}+{P}^{2}\right)+2\gamma Qq\sqrt{1-\frac{{Q}^{2}+{P}^{2}}{4}}-{\omega }_{0}.\end{equation} \tag{ 5 }$$

We define the rescaled energy corresponding to h_cl as

$$\begin{equation}{\epsilon}=E/j,\end{equation} \tag{ 6 }$$

which determines an effective Planck constant ℏ_eff = 1/j [60].

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 ω = ω₀ = 1, coupling parameter in the superradiant phase γ = 2γ_c, and pseudo-spin value j = 30 ($\mathcal{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 [54].

A PO $\mathcal{O}$ with period T is a subset of the phase space, such that

$$\begin{equation}\mathcal{O}=\left\{\boldsymbol{x}\left(t\right)\enspace \vert \enspace t\in \left[0,T\right]\quad \text{and}\quad \boldsymbol{x}\left(0\right)=\boldsymbol{x}\left(T\right)\right\}.\end{equation} \tag{ 7 }$$

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 [59]:

$$\begin{equation*}{\boldsymbol{x}}_{\text{GS}}=\left(q=-\sqrt{\frac{4{\gamma }^{2}}{{\omega }^{2}}-\frac{{\omega }_{0}^{2}}{4{\gamma }^{2}}},Q=\sqrt{2-\frac{\omega {\omega }_{0}}{2{\gamma }^{2}}},p=P=0\right),\end{equation*}$$

and ${\tilde {\boldsymbol{x}}}_{\text{GS}}$, where the signs of q and Q are opposite to those of x _GS. In our case,

$$\begin{equation*}{\boldsymbol{x}}_{\text{GS}}=\left(q=-1.936,Q=1.225,p=P=0\right)\quad \text{and}\quad {\tilde {\boldsymbol{x}}}_{\text{GS}}=\left(q=1.936,Q=-1.225,p=P=0\right).\end{equation*}$$

We first consider x _GS, which is marked with a blue arrow in figure 1 (a1) and (a2), and with a red arrow in figure 1 (b1) and (b2).

Zoom In Zoom Out Reset image size

By considering small displacements around the stationary point x _GS, we obtain two normal frequencies of the system, ${{\Omega}}_{{{\epsilon}}_{\text{GS}}}^{A}$ and ${{\Omega}}_{{{\epsilon}}_{\text{GS}}}^{B}$, which are given by

$$\begin{equation*}{{\Omega}}_{{{\epsilon}}_{\text{GS}}}^{A,B}=\sqrt{\frac{1}{2{\omega }^{2}}\left(\left(16{\gamma }^{4}+{\omega }^{4}\right){\pm}\sqrt{{\left({\omega }^{4}-16{\gamma }^{4}\right)}^{2}+4{\omega }^{6}{\omega }_{0}^{2}}\right)},\end{equation*}$$

where A corresponds to the plus sign and B to the minus sign. For our selection of parameters, ${{\Omega}}_{{{\epsilon}}_{\text{GS}}}^{A}=4.008$ and ${{\Omega}}_{{{\epsilon}}_{\text{GS}}}^{B}=0.966$. Correspondingly, we have two normal periods ${T}_{{{\epsilon}}_{\text{GS}}}^{A}=2\pi /{{\Omega}}_{{{\epsilon}}_{\text{GS}}}^{A}=1.568$ and ${T}_{{{\epsilon}}_{\text{GS}}}^{B}=2\pi /{{\Omega}}_{{{\epsilon}}_{\text{GS}}}^{B}=6.503$.

Let us focus first on the normal mode of period ${T}_{{{\epsilon}}_{\text{GS}}}^{A}$ with the trivial PO ${\mathcal{O}}_{{{\epsilon}}_{\text{GS}}}^{A}=\left\{{\boldsymbol{x}}_{\text{GS}}\right\}$. By perturbing this stable stationary point and using a monodromy method [61] to guarantee the convergence, we can find a new PO ${\mathcal{O}}_{{{\epsilon}}^{\prime }}^{A}$ with energy ' = _GS + δ and period ${T}_{{{\epsilon}}^{\prime }}^{A}={T}_{{{\epsilon}}_{\text{GS}}}^{A}+\delta T$ [62]. This procedure can be successively repeated to increasing energies , so that a continuous family of periodic orbits ${\mathcal{O}}_{{\epsilon}}^{A}$ is obtained all the way to the chaotic regime, where the orbits become unstable (see appendix A for details). This family of POs is denoted by $\mathcal{A}$,

$$\begin{equation}\mathcal{A}=\left\{{\mathcal{O}}_{{\epsilon}}^{A}\enspace \vert \enspace {{\epsilon}}_{\text{GS}}{\leqslant}{\epsilon}{\leqslant}0\right\},\end{equation} \tag{ 8 }$$

and it is shown in figure 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}_{{{\epsilon}}_{\text{GS}}}^{B}$, yielding the POs of family $\mathcal{B}$,

$$\begin{equation}\mathcal{B}=\left\{{\mathcal{O}}_{{\epsilon}}^{B}\enspace \vert \enspace {{\epsilon}}_{\text{GS}}{\leqslant}{\epsilon}{\leqslant}0\right\},\end{equation} \tag{ 9 }$$

which is plotted in figure 1 (b1) and (b2).

The energy of each PO in the two families identified increases as the POs grow away from x _GS. In figure 1 (c), we show the period ${T}_{{\epsilon}}^{A}$ of the POs in $\mathcal{A}$ (blue line) and ${T}_{{\epsilon}}^{B}$ for the POs in $\mathcal{B}$ (red line) as a function of the energy . For a given energy , the period ${T}_{{\epsilon}}^{B}$ is between three and four times larger than ${T}_{{\epsilon}}^{A}$, a feature to which we come back to when explaining our next results. In figure 1 (d), we display the Lyapunov exponents of those POs as a function of energy (see appendix B for details). The orbits become unstable when the Lyapunov exponents are different from zero. Notice that for family $\mathcal{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 $\mathcal{O}$ is a PO, we may find another PO by mirroring

$$\begin{equation*}\begin{gathered}\mathcal{O}\hspace{25.0pt}\quad {\to }\tilde {\mathcal{O}}\quad \\ \left(q,p;Q,P\right){\mapsto}\left(-q,p;-Q,P\right),\end{gathered}\end{equation*}$$

where $\tilde {\mathcal{O}}$ has the same period, energy, and Lyapunov exponent as $\mathcal{O}$. This transformation yields the mirrored families

$$\begin{equation}\tilde {\mathcal{A}}=\left\{{\tilde {\mathcal{O}}}_{{\epsilon}}^{A}\enspace \vert \enspace {\mathcal{O}}_{{\epsilon}}^{A}\in \mathcal{A}\right\}\quad \text{and}\quad \tilde {\mathcal{B}}=\left\{{\tilde {\mathcal{O}}}_{{\epsilon}}^{B}\enspace \vert \enspace {\mathcal{O}}_{{\epsilon}}^{B}\in \mathcal{B}\right\}.\end{equation} \tag{ 10 }$$

The two families $\tilde {\mathcal{A}}$ and $\tilde {\mathcal{B}}$ are the ones that emanate from the stationary point ${\tilde {\boldsymbol{x}}}_{\text{GS}}$.

The four families $\mathcal{A}$, $\mathcal{B}$, $\tilde {\mathcal{A}}$, and $\tilde {\mathcal{B}}$ scar many of the quantum eigenstates, as we show in this subsection.

3.2.1. Measure of scarring degree

The Husimi function of a state $\hat{\rho }$ 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 $\hat{\rho }$,

$$\begin{equation}{\mathcal{Q}}_{\hat{\rho }}\left(\boldsymbol{x}\right)=\left\langle \boldsymbol{x}\right\vert \hat{\rho }\left\vert \boldsymbol{x}\right\rangle ,\end{equation} \tag{ 11 }$$

where $\left\vert \boldsymbol{x}\right\rangle$ is the coherent state centered at x . In the case of a pure state $\hat{\rho }=\left\vert \psi \right\rangle \left\langle \psi \right\vert$, the function reduces to

$$\begin{equation}{\mathcal{Q}}_{\psi }\left(\boldsymbol{x}\right)={\left\vert \left\langle \psi \vert \boldsymbol{x}\right\rangle \right\vert }^{2}.\end{equation} \tag{ 12 }$$

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 $\hat{\rho }$ and the PO $\mathcal{O}$ with period T, we define the quantity

$$\begin{equation}{\left\langle {\mathcal{Q}}_{\hat{\rho }}\right\rangle }_{\mathcal{O}}=\frac{1}{T}{\int }_{0}^{T}\mathrm{d}t\enspace {\mathcal{Q}}_{\hat{\rho }}\left(\boldsymbol{x}\left(t\right)\right),\end{equation} \tag{ 13 }$$

where $\boldsymbol{x}\left(t\right)\in \mathcal{O}$ and any initial point $\boldsymbol{x}\left(0\right)\in \mathcal{O}$ can be used to perform the average. This measure is similar to the tube phase-space projection introduced in reference [13], but now taking into account the temporal behavior of the orbit. From the definition, we have that

$$\begin{equation}{\left\langle {\mathcal{Q}}_{\hat{\rho }}\right\rangle }_{\mathcal{O}}=\text{tr}\left(\hat{\rho }\enspace {\hat{\rho }}_{\mathcal{O}}\right),\end{equation} \tag{ 14 }$$

where

$$\begin{equation}{\hat{\rho }}_{\mathcal{O}}=\frac{1}{T}{\int }_{0}^{T}\mathrm{d}t\enspace \left\vert \boldsymbol{x}\left(t\right)\right\rangle \left\langle \boldsymbol{x}\left(t\right)\right\vert \end{equation} \tag{ 15 }$$

is a tubular Gaussian distribution around the PO and $\left\vert \boldsymbol{x}\left(t\right)\right\rangle$ is the coherent state centered at the point x (t). See appendix C for an illustration of the Husimi function of ${\hat{\rho }}_{\mathcal{O}}$.

From equation (14), we can see that ${\left\langle {\mathcal{Q}}_{\hat{\rho }}\right\rangle }_{\mathcal{O}}$ is the overlap of state $\hat{\rho }$ with ${\hat{\rho }}_{\mathcal{O}}$. To construct a baseline to compare the value of ${\left\langle {\mathcal{Q}}_{\hat{\rho }}\right\rangle }_{\mathcal{O}}$ with, we consider a totally delocalized state

$$\begin{equation}{\hat{\rho }}_{{\epsilon}}=\frac{1}{\mathcal{V}\left({\epsilon}\right)}{\int }_{\mathcal{M}}\mathrm{d}\boldsymbol{x}\enspace \left\vert \boldsymbol{x}\right\rangle \left\langle \boldsymbol{x}\right\vert \delta \left({\epsilon}-{h}_{\text{cl}}\left(\boldsymbol{x}\right)\right)\end{equation} \tag{ 16 }$$

comprised of all the coherent states within a single energy shell at ${\epsilon}={h}_{\text{cl}}\left(\mathcal{O}\right)$, where $\mathcal{V}\left({\epsilon}\right){=\int }_{\mathcal{M}}\hspace{2pt}\mathrm{d}\boldsymbol{x}\delta \left({\epsilon}-{h}_{\text{cl}}\left(\boldsymbol{x}\right)\right)$ is the phase-space volume of the energy shell at . The value of the trace $\text{tr}\left({\hat{\rho }}_{{\epsilon}}\enspace {\hat{\rho }}_{\mathcal{O}}\right)$ gives the overlap between a totally delocalized state and the orbit. By defining

$$\begin{equation}\mathcal{P}\left(\mathcal{O},\hat{\rho }\right)=\frac{\text{tr}\left(\hat{\rho }\enspace {\hat{\rho }}_{\mathcal{O}}\right)}{\text{tr}\left({\hat{\rho }}_{{\epsilon}}{\hat{\rho }}_{\mathcal{O}}\right)},\end{equation} \tag{ 17 }$$

we obtain a direct measure of quantum scarring. A value of $\mathcal{P}\left(\mathcal{O},\hat{\rho }\right)=1$ indicates that the overlap between state $\hat{\rho }$ and the PO $\mathcal{O}$ is equal to that of a totally delocalized state. Values greater than 1 indicate that the state $\hat{\rho }$ is scarred by the PO $\mathcal{O}$. A value of $\mathcal{P}\left(\mathcal{O},\hat{\rho }\right)=2$, for example, says that state $\hat{\rho }$ is twice as likely to be found near the PO as compared to the delocalized state ${\hat{\rho }}_{{\epsilon}}$. Values less than 1 signal that state $\hat{\rho }$ 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 $\hat{\rho }$ is scarred by the POs of the families described in section 3.1 by defining the numbers

$$\begin{equation}{\mathcal{P}}^{A}\left({\epsilon},\hat{\rho }\right)=\mathcal{P}\left({\mathcal{O}}_{{\epsilon}}^{A},\hat{\rho }\right)+\mathcal{P}\left({\tilde {\mathcal{O}}}_{{\epsilon}}^{A},\hat{\rho }\right)\quad \text{and}\quad {\mathcal{P}}^{B}\left({\epsilon},\hat{\rho }\right)=\mathcal{P}\left({\mathcal{O}}_{{\epsilon}}^{B},\hat{\rho }\right)+\mathcal{P}\left({\tilde {\mathcal{O}}}_{{\epsilon}}^{B},\hat{\rho }\right)\end{equation} \tag{ 18 }$$

for any energy . In words, ${\mathcal{P}}^{A}$ measures how much $\hat{\rho }$ populates the region of the phase space visited by the POs of families $\mathcal{A}$ and $\tilde {\mathcal{A}}$ at the energy , while ${\mathcal{P}}^{B}$ does the same for the families $\mathcal{B}$ and $\tilde {\mathcal{B}}$.

3.2.2. Overlap of the energy eigenstates with classical periodic orbits

For an eigenstate ${\hat{\rho }}_{k}=\left\vert {E}_{k}\right\rangle \left\langle {E}_{k}\right\vert$ of the Dicke Hamiltonian ${\hat{H}}_{\mathrm{D}}$ with scaled eigenenergy _k = E_k /j, the numbers

$$\begin{equation}{\mathcal{P}}_{k}^{A}={\mathcal{P}}^{A}\left({{\epsilon}}_{k},{\hat{\rho }}_{k}\right)\quad \text{and}\quad {\mathcal{P}}_{k}^{B}={\mathcal{P}}^{B}\left({{\epsilon}}_{k},{\hat{\rho }}_{k}\right)\end{equation} \tag{ 19 }$$

measure the scarring produced by the orbits of families $\left(\mathcal{A},\tilde {\mathcal{A}}\right)$ and $\left(\mathcal{B},\tilde {\mathcal{B}}\right)$, respectively ⁶ .

In figure 2, we show the expectation value of the operator ${\hat{j}}_{z}={\hat{J}}_{z}/j$ for each energy eigenstate, leading to an arrangement known as a Peres lattice [51, 63]. 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 figure 2(a), we color the points according to the degree of scarring of each eigenstate $\left\vert {E}_{k}\right\rangle$ by the family $\mathcal{A}$, that is, according to the value of ${\mathcal{P}}_{k}^{A}$. The same is done in figure 2(b) for family $\mathcal{B}$ using now ${\mathcal{P}}_{k}^{B}$.

Zoom In Zoom Out Reset image size

As the phase-space distribution of a scarred state is close to a PO, the expectation value of the observable ${\hat{j}}_{z}$ should be similar to the average of the classical variable j_z = (Q² + P²)/2 over the PO. In figures 2(a) and (b), we plot with blue and red thick lines the value of the classical average

$$\begin{equation}{\left\langle {\hspace{1.5pt}j}_{z}\right\rangle }_{{\mathcal{O}}_{{\epsilon}}^{A,B}}=\frac{1}{{T}_{{\epsilon}}^{A,B}}{\int }_{0}^{{T}_{{\epsilon}}^{A,B}}\mathrm{d}t\enspace {j}_{z}\left(\boldsymbol{x}\left(t\right)\right),\end{equation} \tag{ 20 }$$

over each PO from family $\mathcal{A}$ (blue) and $\mathcal{B}$ (red), where ${T}_{{\epsilon}}^{A,B}$ is the period of ${\mathcal{O}}_{{\epsilon}}^{A,B}$. We see in figure 2(a) [figure 2(b)] that for the low energies in the regular regime, the eigenstates with a high value of ${\mathcal{P}}_{k}^{A}$ [${\mathcal{P}}_{k}^{B}$] lie very close to this blue [red] line.

Using the Bohr–Sommerfeld quantization rule, we semiclassically quantize the energies of the families $\mathcal{A}$ and $\mathcal{B}$. Beginning from the ground-state energy ${\mathcal{E}}_{0}^{A}={\mathcal{E}}_{0}^{B}={{\epsilon}}_{\text{GS}}$, we successively find the excited energies ${\mathcal{E}}_{i}^{A,B}$ by applying the quantization condition,

$$\begin{equation}{\int }_{{\mathcal{E}}_{i-1}^{A,B}}^{{\mathcal{E}}_{i}^{A,B}}\mathrm{d}{\epsilon}\enspace {T}_{{\epsilon}}^{A,B}=2\pi {\hslash }_{\text{eff}}=\frac{2\pi }{j}.\end{equation} \tag{ 21 }$$

We plot the obtained energies ${\mathcal{E}}_{i}^{A}$ [${\mathcal{E}}_{i}^{B}$] with vertical thin blue [red] lines in figure 2(a) [figure 2(b)]. These vertical lines coincide perfectly with the individual eigenstates scarred by the POs of the families $\mathcal{A}$ and $\mathcal{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 distributed around the semiclassical energies. According to Gutzwiller trace formula for the density of states [64] and also to reference [3], the width of these clusters should be given by the Lyapunov exponent of the PO at the respective energy, which is indeed confirmed by the horizontal bars shown at the bottom of figures 2(a) and (b), whose width is twice the value of the Lyapunov exponents of ${\mathcal{O}}_{{\mathcal{E}}_{i}^{A}}^{A}$ (blue) and ${\mathcal{O}}_{{\mathcal{E}}_{i}^{B}}^{B}$ (red) multiplied by ℏ_eff = 1/j.

The distribution of scarred states in the spectrum is governed by two numbers: (a) the periods of the POs determine the semiclassically quantized energies ${\mathcal{E}}_{i}$ around which clusters of scarred states appear, and (b) the Lyapunov exponents control the width of these clusters. These two numbers behave differently for families $\mathcal{A}$ and $\mathcal{B}$ as explained below.

• (a)  
  The periods of the POs in family $\mathcal{A}$ are between three and four times shorter than those of $\mathcal{B}$. This has two consequences, ${\mathcal{P}}_{k}^{A}$ is higher than ${\mathcal{P}}_{k}^{B}$ and there is a smaller number of eigenstates scarred by family $\mathcal{A}$ as compared to $\mathcal{B}$. The former occurs because the POs from family $\mathcal{A}$ are shorter in the phase space than those from family $\mathcal{B}$, thus having smaller overlaps with a delocalized state. Since the denominator in $\mathcal{P}\left(\mathcal{O},\hat{\rho }\right)$ in equation (17) is the overlap of the PO with a delocalized state, this results in the values of ${\mathcal{P}}_{k}^{A}$ being between three to four times higher than ${\mathcal{P}}_{k}^{B}$. The latter occurs because a shorter period translates in greater energy separations ${\mathcal{E}}_{i+1}-{\mathcal{E}}_{i}$. Thus, there are approximately four times more red vertical lines in figure 2(b), marking ${\mathcal{E}}_{i}^{B}$, than blue vertical lines in figure 2(a), marking ${\mathcal{E}}_{i}^{A}$. Because scarred eigenstates appear close to these energies, this results in more states scarred by family $\mathcal{B}$ than states scarred by family $\mathcal{A}$.
• (b)  
  The Lyapunov times of the POs in family $\mathcal{A}$ are always much larger than their periods, which reflects the fact that the energy differences ${\mathcal{E}}_{i+1}^{A}-{\mathcal{E}}_{i}^{A}$ are larger than the Lyapunov exponents multiplied by ℏ_eff. For family $\mathcal{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 $\mathcal{A}$ around the semiclassical energies ${\mathcal{E}}_{i}^{A}$ in figure 2(a), but makes it harder to distinguish the ones scarred by $\mathcal{B}$ around ${\mathcal{E}}_{i}^{B}$ in figure 2(b), because neighboring clusters may overlap.

We close section 3.2.2 with the interesting observation that POs from both families may affect the same eigenstate. In figure 2(a) [figure 2(b)], the blue circles outlined by red [red circles outlined by blue] mark the eigenstates where both ${\mathcal{P}}_{k}^{A}$ and ${\mathcal{P}}_{k}^{B}$ are greater than 4. These five states are located at energies where the semiclassical quantizations of both families coincide.

3.2.3. Projected Husimi distribution

We now use the Husimi distribution to visually confirm that the energy eigenstates with large values of ${\mathcal{P}}_{k}^{A}$ and ${\mathcal{P}}_{k}^{B}$ are scarred by the classical POs of families $\mathcal{A}$ and $\mathcal{B}$.

For a state $\hat{\rho }$ and energy , we consider the projection of the Husimi function ${\mathcal{Q}}_{\hat{\rho }}$ over the classical energy shell h_cl( x ) = by integrating out the bosonic variables (q, p),

$$\begin{equation}{\tilde {\mathcal{Q}}}_{{\epsilon},\hat{\rho }}\left(Q,P\right)=\iint \mathrm{d}q\enspace \mathrm{d}p\enspace \delta \left({\epsilon}-{h}_{\text{cl}}\left(\boldsymbol{x}\right)\right){\mathcal{Q}}_{\hat{\rho }}\left(\boldsymbol{x}\right),\end{equation} \tag{ 22 }$$

where x = (q, p; Q, P). For an eigenstate ${\hat{\rho }}_{k}=\left\vert {E}_{k}\right\rangle \left\langle {E}_{k}\right\vert$, the projection ${\tilde {\mathcal{Q}}}_{k}={\tilde {\mathcal{Q}}}_{{{\epsilon}}_{k},{\hat{\rho }}_{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 appendix D for details on the computation of this projection).

We select 12 energy eigenstates [56], marked as A1–A6 and B1–B6 in figures 2(a) and (b), which in addition to having high values of ${\mathcal{P}}_{k}^{A,B}$ lie close to the classical average of j_z given by equation (20). We plot the Husimi projection ${\tilde {\mathcal{Q}}}_{k}\left(Q,P\right)$ for each one of these eigenstates at the 12 bottom panels of figure 2. These distributions are superposed by the projections of the POs ${\mathcal{O}}_{{{\epsilon}}_{k}}^{A}$ (blue solid line), ${\mathcal{O}}_{{{\epsilon}}_{k}}^{B}$ (red solid line), ${\tilde {\mathcal{O}}}_{{{\epsilon}}_{k}}^{A}$ (blue dashed line), and ${\tilde {\mathcal{O}}}_{{{\epsilon}}_{k}}^{B}$ (red dashed line). 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 toward 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 figure 6(b) in appendix C 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 ${\mathcal{P}}_{k}^{A}$ for A5 and A6 as compared to A1–A4.

We expand the selected coherent states in the Hamiltonian eigenbasis

$$\begin{equation}\left\vert \boldsymbol{x}\right\rangle =\sum\limits _{k}{c}_{k}\left\vert {E}_{k}\right\rangle .\end{equation} \tag{ 23 }$$

The amplitudes c_k and the degree of scarring of the corresponding eigenstates $\left\vert {E}_{k}\right\rangle$ determine the dynamical properties of the initial states.

We consider three initial coherent states, $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$, $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$, and $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$, where

$$\begin{equation}\begin{aligned}\hfill {\boldsymbol{x}}_{\text{i}}& =\left(q=-0.795,p=0;Q=1.75,P=0\enspace \right),\hfill \\ \hfill {\boldsymbol{x}}_{\text{ii}}& =\left(q=-0.105,p=0;Q=0.9,P=0.7\right),\hfill \\ \hfill {\boldsymbol{x}}_{\text{iii}}& =\left(q=0.624,p=0;Q=-0.2,P=1\right).\hfill \end{aligned}\end{equation} \tag{ 24 }$$

They all have mean energy in the chaotic region, $\langle {\hat{H}}_{D}\rangle /j={{\epsilon}}_{\text{i}}={{\epsilon}}_{\text{ii}}={{\epsilon}}_{\text{iii}}=-0.5$, and the widths of their energy distributions in the energy eigenbasis are (σ_i, σ_ii, σ_iii) = (0.248, 0.212, 0.192) (in units of ) [32, 65]. The points x _i, x _ii, and x _iii are marked with little circles in figure 1 (a1), (a2), (b1), and (b2). They all have the same shade of green indicating their equal energy. The point x _i sits on top of the PO of family $\mathcal{A}$ that has that same energy −0.5, as can be seen in figure 1 (a1) and (a2). The point x _ii sits on top of the PO of $\mathcal{B}$, as shown in figure 1 (b1) and (b2). The point x _iii is located far away from the POs of both families (and the mirrored families), as seen in figure 1 (a1), (a2), (b1), and (b2).

4.1.1. Scarring of coherent states

In figure 3 (ai)–(aiii) and (bi)–(biii), we plot the energy distribution,

$$\begin{equation}\mathcal{G}\left({\epsilon}\right)=\sum\limits _{k}{\left\vert {c}_{k}\right\vert }^{2}\delta \left({\epsilon}-{{\epsilon}}_{k}\right),\end{equation} \tag{ 25 }$$

of the initial coherent states $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ [figure 3 (ai) and (bi)], $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$ [figure 3 (aii) and (bii)], and $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$ [figure 3 (aiii) and (biii)]. These distributions are usually referred to as local density of states (LDoS). In a sense, figure 3 (ai)–(biii) are similar to figures 2(a) and (b), but now instead of j_z in the vertical axis, we have the components of the chosen coherent states. The distributions in figure 3 (ai)–(aiii) are the same as in figure 3 (bi)–(biii), what changes is just the colors: the blue tones in figure 3 (ai)–(aiii) indicate the values of ${\mathcal{P}}_{k}^{A}$ and the red tones in figure 3 (bi)–(biii) indicate the values of ${\mathcal{P}}_{k}^{B}$ [see equation (19)]. In addition to the circles corresponding to amplitudes of the weights ${\left\vert {c}_{k}\right\vert }^{2}$, figure 3 (ai)–(aiii) [figure 3 (bi)–(biii)] also display vertical lines that mark the semiclassically quantized energies ${\mathcal{E}}_{i}^{A}$ [${\mathcal{E}}_{i}^{B}$] given by equation (21).

Zoom In Zoom Out Reset image size

We observe the following features for the three initial coherent states:

• (a)  
  For the initial coherent state $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ located in the PO of family $\mathcal{A}$, the largest components ${\left\vert {c}_{k}\right\vert }^{2}$, indicated as i1-i6 in figure 3 (ai), correspond to the eigenstates with large values of ${\mathcal{P}}_{k}^{A}$ and therefore scarred by the POs of family $\mathcal{A}$. Due to the high participation of these eigenstates, we say that the initial coherent state $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ is itself scarred by family $\mathcal{A}$. As visible in figure 3 (ai), the LDoS of $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ exhibits a clear comb-like pattern, which is typical of scarred states [3]. As seen in figure 2(a), the scarred eigenstates cluster around the semiclassical energies ${\mathcal{E}}_{i}^{A}$. Because these scarred eigenstates give the largest contributions to the initial state, that is they lead to the biggest components ${\left\vert {c}_{k}\right\vert }^{2}$, the LDoS attains higher values around the energies ${\mathcal{E}}_{i}^{A}$, creating the comb-like structure.One sees that the contributions to $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ from eigenstates non-scarred by the POs of family $\mathcal{A}$ are erratically distributed, with small and medium values of ${\left\vert {c}_{k}\right\vert }^{2}$, as evident from the red points in figure 3 (bi), which mark the eigenstates according to their values of ${\mathcal{P}}_{k}^{B}$.
• (b)  
  A similar picture emerges for the initial coherent state $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$ located in the PO of family $\mathcal{B}$. Its largest components ${\left\vert {c}_{k}\right\vert }^{2}$, indicated as ii1–ii6 in figure 3 (bii), correspond to the eigenstates with large values of ${\mathcal{P}}_{k}^{B}$ and thus scarred by the POs of family $\mathcal{B}$. The contributions from eigenstates non-scarred by this family are smaller and their values fluctuate randomly. The initial coherent state $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$ is therefore scarred by family $\mathcal{B}$. Its comb-like structure is somewhat visible, but, because the separations ${\mathcal{E}}_{i}^{B}-{\mathcal{E}}_{i+1}^{B}$ are of the order of the Lyapunov exponents of the POs in family $\mathcal{B}$, it is harder to distinguish it.As discussed in the previous section, there are more eigenstates scarred by family $\mathcal{B}$ than by family $\mathcal{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 figure 3 (ai) with the number of red circles in figure 3 (bii). The larger number of contributing states from family $\mathcal{B}$ explains why the biggest components in figure 3 (bii) are smaller than the ones in figure 3 (ai).
• (c)  
  For the initial coherent state $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$, which is located far away from the POs of families $\mathcal{A}$, $\mathcal{B}$, $\tilde {\mathcal{A}}$ and $\tilde {\mathcal{B}}$, the largest components ${\left\vert {c}_{k}\right\vert }^{2}$, indicated as iii1–iii6 in figure 3 (aiii) and (biii), are not colored by any of the two families, as expected.

In figure 3 (ci)–(ciii), we plot the same Peres lattice of the pseudo-spin operator ${\hat{j}}_{z}$ that we showed in figure 2, but we now use tones of green to indicate the values of ${\left\vert {c}_{k}\right\vert }^{2}$ for the coherent states $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ [figure 3 (ci)], $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$ [figure 3 (cii)], and $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$ [figure 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 $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ ($\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$) concentrate toward the top (bottom) of the lattice, because this is where the eigenstates scarred by family $\mathcal{A}$ ($\mathcal{B}$) tend to be found, closer to the classical average of j_z (equation (20)) over the corresponding family [blue (red) line in figure 3 (ci)–(ciii)]. The components of $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$ spread more homogeneously over the middle part of the lattice, as seen in figure 3 (ciii).

4.1.2. Husimi distributions of the most contributing eigenstates

In figure 4, we plot the projected Husimi distributions for each of the six eigenstates with the largest components in the LDoS of the initial coherent states $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ [figure 4 (i1)–(i6)], $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$ [figure 4 (ii1)–(ii6)], and $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$ [figure 4 (iii1)–(iii6)]. This is done in two different ways:

Zoom In Zoom Out Reset image size

(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 figure 2 A1–B6 [see equation (22)] and in [56]. This choice is shown in green in the top of each panel of figure 4.

(b) By directly integrating over the bosonic variables (q, p) of the Husimi function, $\iint \mathrm{d}q\enspace \mathrm{d}p\enspace {\mathcal{Q}}_{k}\left(\boldsymbol{x}\right)$, as done in [55, 57, 66]. This alternative is plotted in orange in the bottom of each panel of figure 4.

As expected, the Husimi functions of the eigenstates i1–i6 [ii1–ii6] with the largest components in the coherent state $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ [$\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$] concentrate along the POs of family $\mathcal{A}$ [$\mathcal{B}$] at the corresponding energy. This is more evident with the Husimi functions intersected by the energy shells (top green plots in figure 4) than in the complete projected Husimi functions (bottom orange plots in figure 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 [56] and more details are given in appendix D.

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 $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$ are not smoothly distributed either. As seen in figure 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 $\mathcal{A}$ and $\mathcal{B}$.

To study the dynamics of the three selected initial states, we consider the survival probability defined as

$$\begin{equation}{S}_{\mathrm{P}}\left(t\right)=\vert \langle \boldsymbol{x}\vert \hat{U}\left(t\right)\vert \boldsymbol{x}\rangle {\vert }^{2},\end{equation} \tag{ 26 }$$

where $\hat{U}\left(t\right)={\text{e}}^{-\text{i}\hspace{1pt}{\hat{H}}_{\mathrm{D}}t}$ is the unitary evolution operator. S_P(t) measures the probability of the evolved state $\hat{U}\left(t\right)\left\vert \boldsymbol{x}\right\rangle$ to return to its initial state $\left\vert \boldsymbol{x}\right\rangle$ for any given time. By expanding the coherent state in the Hamiltonian eigenbasis as in equation (23), we can rewrite the survival probability as

$$\begin{equation}{S}_{\mathrm{P}}\left(t\right)={\left\vert \sum\limits _{k}{\left\vert {c}_{k}\right\vert }^{2}\enspace {\mathrm{e}}^{-\mathrm{i}{{\epsilon}}_{k}t/{\hslash }_{\text{eff}}}\right\vert }^{2},\end{equation} \tag{ 27 }$$

which can also be written in an integral form,

$$\begin{equation}{S}_{\mathrm{P}}\left(t\right)={\left\vert \int \mathrm{d}{\epsilon}\enspace \mathcal{G}\left({\epsilon}\right){\text{e}}^{-\text{i}{\epsilon}t/{\hslash }_{\text{eff}}}\right\vert }^{2}={\left\vert \mathcal{F}\left[\mathcal{G}\right]\left(t\right)\right\vert }^{2},\end{equation} \tag{ 28 }$$

which is the squared norm of the Fourier transform of the LDoS, as defined in equation (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 values of the survival probability are larger than those for non-scarred states [13, 36]. 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 $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$, $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$, and $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ are respectively shown in figure 5 (ai)–(bi), figure 5 (aii)–(bii), and figure 5 (aiii)–(biii). The orange curves are numerical results and the thick purple 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 [36]. The subsequent behavior differs among the three states.

Zoom In Zoom Out Reset image size

For state $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$, the first revival of the survival probability occurs at t = 2.8 [figure 5 (bi)], which is precisely the average of the periods ${T}_{{\epsilon}}^{A}$ within the energy range of state $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$, ${T}_{\text{i}}^{A}$. A second revival is observed at $t=2{T}_{\text{i}}^{A}$, while the next ones do not follow this period, because ${T}_{{\epsilon}}^{A}$ changes strongly through the energy range of state $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ destroying the periodicity. This may be seen in the left inset of figure 1 (c), where the little blue square marks the value of ${T}_{\text{i}}^{A}$. For longer times, the survival probability equilibrates, fluctuating around its asymptotic value drawn with a horizontal black dashed line [36].

The first revival of S_P(t) for state $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$ appears at $t={T}_{\text{ii}}^{B}=9.5$ [figure 5 (bii)], which corresponds to the average period of the POs in the family $\mathcal{B}$ within the energy range of state $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$. The subsequent revivals happens at multiples of this period. The periods ${T}_{{\epsilon}}^{B}$ do not vary as much as ${T}_{{\epsilon}}^{A}$ around = −0.5, as seen in the right inset of figure 1 (c), where the value ${T}_{\text{ii}}^{B}$ is marked with a little red square. The slope of ${T}_{{\epsilon}}^{A}$ around ${T}_{\text{i}}^{A}$ is large, while ${T}_{{\epsilon}}^{B}$ actually attains a local minimum close to ${T}_{\text{ii}}^{B}$. As a result, the revivals of the survival probability of state $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$ follow the period ${T}_{\text{ii}}^{B}$ for much longer than in the case of state $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$.

A great advantage of having identified the families $\mathcal{A}$ and $\mathcal{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 $\mathcal{B}$ than from family $\mathcal{A}$ help us understand also why the survival probability for the initial state $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$ saturates at a lower point than S_P(t) for the state $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$.

In contrast to states $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ and $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$, the initial coherent state $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$ 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 figure 5 (aiii) and (biii) and it is well described using random matrix theory [33, 36].

A complete picture of the frequencies dominating the dynamical behavior of the coherent states may be obtained by calculating the inverse Fourier transform of the survival probability,

$$\begin{equation}{\mathcal{F}}^{-1}\left[{S}_{\mathrm{P}}\right]\left({\epsilon}\right)=\frac{1}{2\pi {\hslash }_{\text{eff}}}\int \mathrm{d}t\enspace {S}_{\mathrm{P}}\left(t\right){\text{e}}^{\text{i}{\epsilon}t/{\hslash }_{\text{eff}}}=\int \mathrm{d}{{\epsilon}}^{\prime }\mathcal{G}\left({{\epsilon}}^{\prime }\right)\enspace \mathcal{G}\left({{\epsilon}}^{\prime }+{\epsilon}\right),\end{equation} \tag{ 29 }$$

which is the autocorrelation function of the LDoS. This function gives us the distribution of frequencies of the survival probability in units of . In figure 5 (ci), (cii), and (ciii), we plot ${\mathcal{F}}^{-1}\left[{S}_{\mathrm{P}}\right]\left({\epsilon}\right)$ with black bars for states $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$, $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$, and $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$, respectively. The vertical blue [red] lines in figure 5 (ci) [figure 5 (cii)] mark the integer multiples of the energy separation $2\pi {\hslash }_{\text{eff}}/{T}_{\text{i}}^{A}$ [$2\pi {\hslash }_{\text{eff}}/{T}_{\text{ii}}^{B}$] corresponding to the period ${T}_{\text{i}}^{A}$ [${T}_{\text{ii}}^{B}$]. In figure 5 (ci) [figure 5 (cii)], the function ${\mathcal{F}}^{-1}\left[{S}_{\mathrm{P}}\right]\left({\epsilon}\right)$ displays a comb-like pattern following the energy separation above, indicating that the period ${T}_{\text{i}}^{A}$ [${T}_{\text{ii}}^{B}$] indeed dominates the behavior of the survival probability of state $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ [$\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$]. Figure 5 (ciii), on the other hand, does not exhibit any recognizable pattern. Because ${\mathcal{F}}^{-1}\left[{S}_{\mathrm{P}}\right]\left({\epsilon}\right)$ is the autocorrelation function of the LDoS (equation (29)), it enhances any periodic structure present in $\mathcal{G}\left({\epsilon}\right)$. This amplifies the comb-like structures that is present in the LDoS of states $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ and $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$, and confirms that such a structure is absent in the LDoS of state $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$.

The fact that the survival probability of state $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$ is so well described by random matrix theory and that no frequencies clearly stand out in its Fourier transform is puzzling at first sight, since figure 4 (iii1)–(iii6) suggest that the most contributing eigenstates to $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$, that is eigenstates iii1–iii6 in figure 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.

We finally analyze how the scarring manifests in the infinite-time average of the coherent states [56]

$$\begin{equation}{\bar{\rho }}_{\boldsymbol{x}}={\mathrm{lim}}_{T\to \infty }\enspace \frac{1}{T}{\int }_{0}^{T}\mathrm{d}t\enspace \hat{U}\left(t\right)\left\vert \boldsymbol{x}\right\rangle \left\langle \boldsymbol{x}\right\vert {\hat{U}}^{{\dagger}}\left(t\right).\end{equation} \tag{ 30 }$$

With the functions ${\mathcal{P}}^{A}$ and ${\mathcal{P}}^{B}$ defined by equation (18), we may calculate

$$\begin{align*}\hfill {\mathcal{P}}_{\text{i}}^{A}& ={\mathcal{P}}^{A}\left({{\epsilon}}_{\text{i}},{\bar{\rho }}_{{\boldsymbol{x}}_{\text{i}}}\right)=4.77\hfill & \hfill {\mathcal{P}}_{\text{ii}}^{A}& ={\mathcal{P}}^{A}\left({{\epsilon}}_{\text{ii}},{\bar{\rho }}_{{\boldsymbol{x}}_{\text{ii}}}\right)=0.96\hfill & \hfill \quad {\mathcal{P}}_{\text{iii}}^{A}& ={\mathcal{P}}^{A}\left({{\epsilon}}_{\text{iii}},{\bar{\rho }}_{{\boldsymbol{x}}_{\text{iii}}}\right)=1.10\hfill & \hfill \\ \hfill {\mathcal{P}}_{\text{i}}^{B}& ={\mathcal{P}}^{B}\left({{\epsilon}}_{\text{i}},{\bar{\rho }}_{{\boldsymbol{x}}_{\text{i}}}\right)=0.78\hfill & \hfill {\mathcal{P}}_{\text{ii}}^{B}& ={\mathcal{P}}^{B}\left({{\epsilon}}_{\text{ii}},{\bar{\rho }}_{{\boldsymbol{x}}_{\text{ii}}}\right)=2.46\hfill & \hfill {\mathcal{P}}_{\text{iii}}^{B}& ={\mathcal{P}}^{B}\left({{\epsilon}}_{\text{iii}},{\bar{\rho }}_{{\boldsymbol{x}}_{\text{iii}}}\right)=0.86.\hfill & \hfill \end{align*}$$

Notice that ${\mathcal{P}}_{\text{i}}^{A}=4.77$ and ${\mathcal{P}}_{\text{ii}}^{B}=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 $\hat{U}\left(t\right)\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ [$\hat{U}\left(t\right)\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$] is more likely to be found in the vicinity of the orbits of families $\mathcal{A}$ and $\tilde {\mathcal{A}}$ [$\mathcal{B}$ and $\tilde {\mathcal{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 [58].

We can visualize the dynamical scars by calculating the projected Husimi distributions ${\mathcal{Q}}_{{\epsilon},\hat{\rho }}$ [equation (22)] of the infinite-time averages $\hat{\rho }={\bar{\rho }}_{{\boldsymbol{x}}_{\text{i}}}$, ${\bar{\rho }}_{{\boldsymbol{x}}_{\text{ii}}}$, and ${\bar{\rho }}_{{\boldsymbol{x}}_{\text{iii}}}$. These are plotted in figure 5 (di), (dii), and (diii), respectively. The dark concentrations in figure 5 (di) follow the orbit of family $\mathcal{A}$ ($\tilde {\mathcal{A}}$) at energy = −0.5 shown with a blue solid (dashed) line. The concentrations in figure 5 (dii) follow the orbit from $\mathcal{B}$ ($\tilde {\mathcal{B}}$) at the same energy shown with a solid (dashed) red line. Interestingly, some dark concentrations are visible also in figure 5 (diii), which suggest that state $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$ is scarred, but by POs from families other than those identified in this work.

The fact that the survival probability of state $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$ is so well described by random matrix theory [green line in figure 5 (aiii) and (biii)] and no identifiable structure is visible in the autocorrelation of its LDoS [figure 5 (ciii) and equation (29))] suggests that even though this state has a minor degree of scarring that is visible in phase space [figure 5 (diii)], it is not identifiable by just looking in the Hilbert space. To better understand this feature, we use the ratio between the asymptotic value of the survival probability ${S}_{\mathrm{P}}^{\infty }$ of $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$ and the asymtpotic value of the survival probability of an ensemble of random states with the same enveloping LDoS ${S}_{\mathrm{P}}^{\left(r\right),\infty }$ [36],

$$\begin{equation}R=\frac{{S}_{\mathrm{P}}^{\left(r\right),\infty }}{{S}_{\mathrm{P}}^{\infty }}.\end{equation} \tag{ 31 }$$

The ratio R measures the similarity between the LDoS of the coherent state and that of a random state, (see reference [36] 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 [33]. For states $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ and $\left\vert {\boldsymbol{x}}_{\text{ii}}\right\rangle$, R_i = 0.41 and R_ii = 0.81, respectively, while R_iii = 1.01 for state $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$. These numbers show that the scarring can be traced back to a comb-like structure in the LDoS of states $\left\vert {\boldsymbol{x}}_{\text{i}}\right\rangle$ and $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$, but not of state $\left\vert {\boldsymbol{x}}_{\text{iii}}\right\rangle$. This is consistent with reference [56], 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 figure 5 (diii). 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.

Assume we have guesses $\check {\boldsymbol{x}}$ and $\check {T}$ for an initial condition and period of a PO, respectively. We want $\boldsymbol{x}=\check {\boldsymbol{x}}+{\Delta}\boldsymbol{x}$ and $T=\check {T}+{\Delta}T$, the initial condition and period of an orbit in the same energy shell as $\check {\boldsymbol{x}}$. 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 [69] 1.1.3). If ${\Vert}{\Delta}\boldsymbol{x}{\Vert}$ is small, one may approximate to first order,

$$\begin{equation*}{\boldsymbol{\varphi }}^{\check {T}}\left(\check {\boldsymbol{x}}+{\Delta}\boldsymbol{x}\right)\approx {\boldsymbol{\varphi }}^{\check {T}}\left(\check {\boldsymbol{x}}\right)+{\mathbf{\Phi }}_{\check {\boldsymbol{x}}}\left(\check {T}\right){\Delta}\boldsymbol{x}.\end{equation*}$$

Similarly, if $\left\vert {\Delta}T\right\vert$ is small, Taylor expanding the flow up to first order,

$$\begin{equation*}{\boldsymbol{\varphi }}^{{\Delta}T}\left(\boldsymbol{x}\right)\approx \boldsymbol{x}+{\Delta}T\enspace \mathbf{\Sigma }\nabla {h}_{\text{cl}}\left(\boldsymbol{x}\right),\quad \mathbf{\Sigma }=\left(\begin{matrix}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \end{matrix}\right),\end{equation*}$$

where ∇h_cl( x ) is the gradient of the Hamiltonian. Then, we have

$$\begin{align}\hfill {\boldsymbol{\varphi }}^{T}\left(\boldsymbol{x}\right)& ={\boldsymbol{\varphi }}^{{\Delta}T+\check {T}}\left(\boldsymbol{x}\right)={\boldsymbol{\varphi }}^{{\Delta}T}\left({\boldsymbol{\varphi }}^{\check {T}}\left(\check {\boldsymbol{x}}+{\Delta}\boldsymbol{x}\right)\right)\hfill \\ \hfill & \approx {\boldsymbol{\varphi }}^{\check {T}}\left(\check {\boldsymbol{x}}+{\Delta}\boldsymbol{x}\right)+{\Delta}T\enspace \mathbf{\Sigma }\nabla {h}_{\text{cl}}\left({\boldsymbol{\varphi }}^{\check {T}}\left(\check {\boldsymbol{x}}+{\Delta}\boldsymbol{x}\right)\right)\hfill \\ \hfill & \approx {\boldsymbol{\varphi }}^{\check {T}}\left(\check {\boldsymbol{x}}\right)+{\mathbf{\Phi }}_{\check {\boldsymbol{x}}}\left(\check {T}\right){\Delta}\boldsymbol{x}+{\Delta}T\enspace \mathbf{\Sigma }\nabla {h}_{\text{cl}}\left({\boldsymbol{\varphi }}^{\check {T}}\left(\check {\boldsymbol{x}}\right)+{\mathbf{\Phi }}_{\check {\boldsymbol{x}}}\left(\check {T}\right){\Delta}\boldsymbol{x}\right)\hfill \\ \hfill & \approx {\boldsymbol{\varphi }}^{\check {T}}\left(\check {\boldsymbol{x}}\right)+{\mathbf{\Phi }}_{\check {\boldsymbol{x}}}\left(\check {T}\right){\Delta}\boldsymbol{x}+{\Delta}T\enspace \mathbf{\Sigma }\nabla {h}_{\text{cl}}\left({\check {\boldsymbol{x}}}^{\prime }\right)\hfill \end{align} \tag{ A1 }$$

where ${\check {\boldsymbol{x}}}^{\prime }={\boldsymbol{\varphi }}^{\check {T}}\left(\check {\boldsymbol{x}}\right)$. Thus, we can approximate the periodicity constriction x = φ ^T ( x ) to first order by

$$\begin{equation}\check {\boldsymbol{x}}+{\Delta}\boldsymbol{x}={\boldsymbol{\varphi }}^{\check {T}}\left(\check {\boldsymbol{x}}\right)+{\mathbf{\Phi }}_{\check {\boldsymbol{x}}}\left(\check {T}\right){\Delta}\boldsymbol{x}+{\Delta}T\enspace \mathbf{\Sigma }\nabla {h}_{\text{cl}}\left({\check {\boldsymbol{x}}}^{\prime }\right).\end{equation} \tag{ A2 }$$

Also, we can approximate the energy constriction ${h}_{\text{cl}}\left(\boldsymbol{x}\right)={h}_{\text{cl}}\left(\check {\boldsymbol{x}}\right)$ to first order to get

$$\begin{equation}\nabla {h}_{\text{cl}}\left({\check {\boldsymbol{x}}}^{\prime }\right)\cdot {\Delta}\boldsymbol{x}=0,\end{equation} \tag{ A3 }$$

and, finally, the constriction to stay in the same Poincaré section of constant P as

$$\begin{equation}\boldsymbol{\xi }\cdot {\Delta}\boldsymbol{x}=0,\end{equation} \tag{ A4 }$$

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

$$\begin{equation*}\left(\begin{matrix}\hfill \left(\boldsymbol{1}-{\mathbf{\Phi }}_{\check {\boldsymbol{x}}}\left(\check {T}\right)\right)\hfill & \hfill -\mathbf{\Sigma }\nabla {h}_{\text{cl}}\left({\check {\boldsymbol{x}}}^{\prime }\right)\hfill \\ \hfill \nabla {h}_{\text{cl}}{\left({\check {\boldsymbol{x}}}^{\prime }\right)}^{\top }\hfill & \hfill 0\hfill \\ \hfill {\boldsymbol{\xi }}^{\top }\hfill & \hfill 0\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {\Delta}\boldsymbol{x}\hfill \\ \hfill {\Delta}T\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {\check {\boldsymbol{x}}}^{\prime }-\check {\boldsymbol{x}}\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{matrix}\right).\end{equation*}$$

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 guesses $\check {\boldsymbol{x}}+{\Delta}\boldsymbol{x}$ and $\check {T}+{\Delta}T$ that will converge to x and T if the initial guesses were good enough.

We now define an iterative process to find the families of POs. For the first step, we start with a stable stationary point ${\boldsymbol{x}}_{\text{GS}}\in \mathcal{M}$ with energy _GS = h_cl( x _GS) and normal period ${T}_{{{\epsilon}}_{\text{GS}}}$.

Given ${\mathcal{O}}_{{\epsilon}}$, we now explain how to find ${\mathcal{O}}_{{{\epsilon}}^{\prime }}$ with ' = δ + close to . Pick $\boldsymbol{x}\in \mathcal{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. Set

$$\begin{equation}{\check {\boldsymbol{x}}}^{\prime }=\boldsymbol{x}+\delta \boldsymbol{x}\quad \text{and}\quad {\check {T}}^{\prime }=T+\delta T,\end{equation} \tag{ A5 }$$

where

$$\begin{equation*}\delta {T}_{k}=\begin{cases}0\quad \hfill & \;\text{if}\;\text{we}\;\text{are}\;\text{in}\;\text{the}\;\text{first}\;\text{step,}\hfill \\ \left({{\epsilon}}^{\prime }-{\epsilon}\right)\frac{T-{T}_{\text{prev}}}{{\epsilon}-{{\epsilon}}_{\text{prev}}}\quad \hfill & \text{else,}\hfill \end{cases}\end{equation*}$$

where _prev and T_prev are the energy and period of the closest previously calculated orbit. This way, ${\check {T}}^{\prime }$ is linear extrapolation based on the behavior of the previous orbits.

Using the monodromy method detailed in the previous subsection, we may correct the guesses ${\check {\boldsymbol{x}}}^{\prime }$ and ${\check {T}}^{\prime }$ to obtain actual solutions x ' and T' so that the desired PO is

$$\begin{equation}{\mathcal{O}}_{{{\epsilon}}^{\prime }}=\left\{{\boldsymbol{\varphi }}^{t}\left({\boldsymbol{x}}^{\prime }\right)\enspace \vert \enspace t\in \left[0,{T}^{\prime }\right]\right\}.\end{equation} \tag{ A6 }$$

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 $\left\vert \delta {\epsilon}\right\vert$ which converge to the desired energy.