Dynamical signatures of quantum chaos and relaxation timescales in a spin-boson system

TL;DR

This paper investigates the correlation hole in the survival probability of the chaotic spin-boson (Dicke) model, deriving an analytical expression and revealing that relaxation times scale linearly with system size, with implications for quantum chaos.

Contribution

It provides an analytical description of the survival probability's evolution in the Dicke model and identifies universal dynamics beyond the correlation hole.

Findings

Analytical expression matches numerical results.

Relaxation time scales linearly with system size.

Dynamics beyond the correlation hole are universal.

Abstract

Quantum systems whose classical counterparts are chaotic typically have highly correlated eigenvalues and level statistics that coincide with those from ensembles of full random matrices. A dynamical manifestation of these correlations comes in the form of the so-called correlation hole, which is a dip below the saturation point of the survival probability's time evolution. In this work, we study the correlation hole in the spin-boson (Dicke) model, which presents a chaotic regime and can be realized in experiments with ultracold atoms and ion traps. We derive an analytical expression that describes the entire evolution of the survival probability and allows us to determine the timescales of its relaxation to equilibrium. This expression shows remarkable agreement with our numerical results. While the initial decay and the time to reach the minimum of the correlation hole depend on the initial state, the dynamics beyond the hole up to equilibration is universal. We find that the relaxation time of the survival probability for the Dicke model increases linearly with system size.

Figures

Equations

$${\hat{H}}_D=\omega {\hat{a}}^{†}\hat{a}+{\omega }_0{\hat{J}}_z+\frac{2\gamma }{\sqrt{\mathcal{N}}}{\hat{J}}_x({\hat{a}}^{†}+\hat{a}).$$

$$\nu(E)=\frac{j}{\omega}\left\{\begin{array}[]{l}\frac{1}{\pi}\int_{y_{-}}^{y_{+}}dy\cos^{-1}\left(\sqrt{\frac{2(y-\epsilon)}{\bar{\gamma}^{2}(1-y^{2})}}\right),\epsilon_{gs}\leq\epsilon<-1\\ \frac{1+\epsilon}{2}+\frac{1}{\pi}\int_{\epsilon}^{y_{+}}dy\cos^{-1}\left(\sqrt{\frac{2(y-\epsilon)}{\bar{\gamma}^{2}(1-y^{2})}}\right),|\epsilon|\leq 1\\ 1,\epsilon>1\end{array}\right.$$

$$S_P(t)=|⟨\Psi (0)|\Psi (t)⟩{|}^2.$$

$$S_P(t)={\left|\sum _k|c_k{|}^2{e}^{-i{E}_{k}t}\right|}^2.$$

$$|c_k{|}^2=\frac{r_{k}f(E_k)}{{\sum }_q{r}_{q}f(E_q)}.$$

$${\rho }^R(E)$$

$${\rho }^G(E)$$

$$\mathcal{C}=\frac{1}{2}\left[\text{erf}\left(\frac{E_{\text{c}}-E_{\text{min}}}{\sqrt{2}{\sigma }_G}\right)-\text{erf}\left(\frac{E_{\text{c}}-E_{\text{max}}}{\sqrt{2}{\sigma }_G}\right)\right],$$

$$S_P(t)={\left|\int {\rho }_0(E)e^{-iEt}dE\right|}^2,$$

$$S_P^R(t)=\frac{{\sin }^2({\sigma }_{R}t)}{({\sigma }_{R}t{)}^2},$$

$$S_P^G(t)=\frac{e^{-{\sigma }_G^2{t}^2}}{{4\mathcal{C}}^2}\mathcal{F}(t),\text{with }$$

$$\mathcal{F}(t)={\left|\text{erf}\left(\frac{E_{\text{c}}-E_{\text{min}}-i{\sigma }_G^{2}t}{\sqrt{2}{\sigma }_G}\right)-\text{erf}\left(\frac{E_{\text{c}}-E_{\text{max}}-i{\sigma }_G^{2}t}{\sqrt{2}{\sigma }_G}\right)\right|}^2.$$

$$S_P^G(t\gg {\sigma }_G^{-1})\approx \frac{1}{2\pi {\mathcal{C}}^2{\sigma }_G^2{t}^2}\times$$

$$\left\{\mathcal{E}-2e^{-\frac{\left[(E_c-E_{\text{min}}{)}^2+(E_{\text{max}}-E_c{)}^2\right]}{2{\sigma }_G^2}}\cos [(E_{\text{max}}-E_{\text{min}})t]\right\},$$

$$\mathcal{E}=\exp \left[-\frac{(E_{\text{c}}-E_{\text{min}}{)}^2}{{\sigma }_G^2}\right]+\exp \left[-\frac{(E_{\text{c}}-E_{\text{max}}{)}^2}{{\sigma }_G^2}\right].$$

$${\bar{S}}_P=\underset{t\to \infty }{\lim }\frac{1}{t}{\int }_0^{t}d{t}^{\prime }{S}_P(t^{\prime }),$$

$$S_P(t)=\sum _{k\ne l}|c_l{|}^2|c_k{|}^2{e}^{-i(E_k-E_l)t}+I_{PR},$$

$$I_{PR}\equiv \sum _k|⟨{\varphi }_k|\Psi (0)⟩{|}^4={\bar{S}}_P$$

$$⟨I_{PR}⟩=⟨\frac{{\sum }_k{r}_k^2{f}^2(E_k)}{{\left({\sum }_q{r}_{q}f(E_q)\right)}^2}⟩\approx \frac{⟨r_k^2⟩}{⟨r_q{⟩}^2}\frac{1}{{\nu }_c}\int {\rho }^2(E)dE,$$

$$⟨I_{PR}^R⟩=\frac{2}{3{\sigma }_R{\nu }_c},$$

$$⟨I_{PR}^G⟩=\frac{\text{erf}\left(\frac{E_{\text{c}}-E_{\text{min}}}{{\sigma }_G}\right)-\text{erf}\left(\frac{E_{\text{c}}-E_{\text{max}}}{{\sigma }_G}\right)}{3\sqrt{\pi }{\sigma }_G{\nu }_c{\mathcal{C}}^2}.$$

$$⟨S_P(t)⟩=\frac{1-⟨I_{PR}⟩}{\eta -1}\left[\eta S_P^{bc}(t)-b_2\left(\frac{t}{2\pi {\nu }_c}\right)\right]+⟨I_{PR}⟩.$$

$$\eta \equiv \frac{{\nu }_c}{\int {\rho }^2(E)dE}=\frac{⟨r_k^2⟩}{⟨r_k{⟩}^2}\frac{1}{⟨I_{PR}⟩},$$

$$b_2(\bar{t})$$

$${\frac{d⟨S_P(t)⟩}{dt}|}_{t=t_{Th}}=0.$$

$$\eta S_P^{bc,R}(t\gg {\sigma }_R^{-1})\to \frac{\eta }{2{\sigma }_R^2{t}^2}.$$

$$e^{-(E_{max}-E_c{)}^2/{\sigma }^2}>\frac{1}{\eta }and{e}^{-(E_{min}-E_c{)}^2/{\sigma }^2}>\frac{1}{\eta }$$

$$\eta S_P^{bc,G}(t\gg {\sigma }_G^{-1})\to \frac{\eta \mathcal{E}}{2\pi {\mathcal{C}}^2{\sigma }_G^2{t}^2}.$$

$$\eta S_P^{bc,G}(t\gg {\sigma }_G^{-1})\to \eta \exp [-{\sigma }_G^2{t}^2].$$

$$b_2\left(\frac{t}{2\pi {\nu }_c}\right)\to 1-\frac{t}{\pi {\nu }_c}\text{for}\frac{t}{{\nu }_c}\ll 1.$$

$$t_{Th}^R={\left(\frac{2\pi {\nu }_c^2}{{\sigma }_R}\right)}^{1/3},$$

$$t_{Th}^G={\left(\frac{\eta {\nu }_c\mathcal{E}}{{\mathcal{C}}^2{\sigma }_G^2}\right)}^{1/3}$$

Full text

I Introduction

The subject of equilibration and thermalization of isolated quantum systems in the chaotic regime has seen a great deal of advance in the last years Reimann (2008); Short (2011); Short and Farrelly (2012); Zangara et al. (2013); He et al. (2013); Gogolin and Eisert (2016); Borgonovi et al. (2016); Alessio et al. (2016); Dymarsky (a); Reimann (2018a); *Reimann2018b. Equilibration is reached when, after a transient time, the observable under investigation shows only small fluctuations around an asymptotic value, and these fluctuations decrease with system size. Thermalization implies that this infinite-time average is very close to the predictions from statistical mechanics, and the difference between the two also decreases with system size. In this picture, an important open question is how long it takes for isolated quantum systems to reach equilibrium. Despite the increasing number of recent works addressing this issue Monnai (2013); Goldstein et al. (2013); Malabarba et al. (2014); Goldstein et al. (2015); Gogolin and Eisert (2016); Reimann (2016); García-Pintos et al. (2017); de Oliveira et al. (2018); Dymarsky (b); Chan et al. (2018); Bertini et al. (2018); Schiulaz et al. ; Borgonovi et al. (2019), there is no agreement regarding how the relaxation timescale should depend on system size, range of interactions, observables, and initial states.

The last steps of the evolution toward equilibrium, after the dynamics resolves the discreteness of the spectrum, are determined by the properties of the eigenvalues Schiulaz et al. . The largest possible timescale is the Heisenberg time Berry (1983), which is proportional to the inverse of the mean level spacing of the region of the spectrum probed by the initial state. Before reaching this timescale, effects of the correlations between the eigenvalues may be observed. In the case of the survival probability footSurv , which is the probability of finding the initial state at time $t$, these correlations cause a decay below the saturation value of the dynamics, known as correlation hole footNote . This phenomenon was first studied in the context of molecules, where the interest was not exactly in dynamics, but in alternative ways to detect level repulsion in systems without good line resolution Leviandier et al. (1986); Guhr and Weidenmüller (1990); Wilkie and Brumer (1991).

The correlation hole has been studied in full random matrices Alhassid and Levine (1992), in many-body systems with Torres-Herrera and Santos (2017a, b); Torres-Herrera et al. (2018); Torres-Herrera and Santos (2019); Schiulaz et al. and without disorder Torres-Herrera and Santos (2017b), in the Sachdev-Ye-Kitaev model Cotler et al. (2017); Gharibyan et al. (2018); Nosaka et al. (2018), which is a two-body random ensemble Brody et al. (1981), and in the finite one-dimensional Anderson model TorresARXIV . The hole is not exclusive to the survival probability, but emerges also in experimental local observables Torres-Herrera et al. (2018); Torres-Herrera and Santos (2019). For the correlation hole to be visible, one needs to perform large averages over initial states and, in the case of Hamiltonian matrices with random elements, over ensembles of Hamiltonian realizations. In Ref. Schiulaz et al. , it was shown that in realistic chaotic many-body quantum systems with local short-range interactions and perturbed far from equilibrium, the time to reach the minimum of the correlation hole increases exponentially with system size. This timescale, which is still shorter than the Heisenberg time, was referred to as Thouless time due to its relationship with the Thouless energy computed from random matrix theory. As explained in Schiulaz et al. , the Thouless time in interacting systems is the time that it takes for an initial state to spread over the entire Hilbert space accessible to its energy. Beyond this point, the dynamics becomes universal all the way to equilibrium.

In the present work, we use the survival probability to study the correlation hole in the Dicke model. This is a paradigmatic spin-boson model with two degrees of freedom. It has a classical counterpart and exhibits chaos for several values of its parameters, mostly for high excitation energies and in the superradiant phase. The model was first introduced to explain superradiance Dicke (1954); Garraway (2011); Baumann et al. (2010, 2011) and has since then been used in different contexts, from quantum chaos Lewenkopf et al. (1991); Emary and Brandes (2003a); *Emary2003; Bastarrachea-Magnani et al. (2014a); *Bastarrachea2015; *Bastarrachea2016PRE; Chávez-Carlos et al. (2016) and quantum batteries Andolina et al. to excited-state quantum phase transitions Pérez-Fernández et al. (2011a); Brandes (2013); Bastarrachea-Magnani et al. (2014b); Larson and Irish (2017) and quench dynamics Pérez-Fernández et al. (2011b); Altland and Haake (2012); Lerma-Hernández et al. (2018); Kloc et al. (2018). Recently, the model was employed in a study of the out-of-time ordered correlator (OTOC), where it was shown that, in the chaotic regime, the OTOC increases exponentially in time with a rate comparable to the classical Lyapunov exponent Chávez-Carlos et al. (2019). In addition to ultracold atoms in optical cavities Baden et al. (2014); Klinder et al. (2015), the Dicke model can now be realized also with ion traps Cohn et al. (2018). The latter is one of the main platforms to study long-time coherence evolution Blatt and Roos (2012), which makes the analysis of the timescales involved in the relaxation process of the Dicke model a timely subject.

We obtain an analytical expression that describes the entire evolution of the survival probability for the Dicke model in the chaotic regime. The expression describes very accurately our numerical results, and with it, we can derive analytically the timescales involved in the relaxation process. We find that the relaxation time increases linearly with system size, while the Thouless time depends non trivially on the initial state.

The article is organized as follows. In Sec. II, we describe the Dicke model and the properties associated with its eigenvalues. In Sec. III, we discuss the initial states considered and present the analytical expression for the survival probability. This expression is compared with numerical results in Sec. IV. The analytical expressions for the Thouless and relaxation times are given and discussed in Sec. V. We present our conclusions in Sec. VI.

II Dicke Model

The Dicke model Dicke (1954) describes the interaction between a set of $\mathcal{N}$ two-level atoms with energy splitting $\omega_{0}$ and a single mode of the electromagnetic field with radiation frequency $\omega$. By setting $\hbar=1$, the time-reversal symmetric Hamiltonian of the model is written as

[TABLE]

The first term of the equation above accounts for the energy of the field, where $\hat{a}^{\dagger}$ ($\hat{a}$) is the bosonic creation (annihilation) operator. The second term corresponds to the energy of the atoms, where $\hat{J}_{x,y,z}=\frac{1}{2}\sum_{k=1}^{\mathcal{N}}\sigma_{x,y,z}^{k}$ are the atomic pseudo-spin operators and $\sigma_{x,y,z}$ are the Pauli matrices. The third term describes the atom-field interaction with coupling parameter $\gamma$. The eigenvalues $j(j+1)$ of the operator $\hat{\textbf{J}}^{2}=\hat{J}_{x}^{2}+\hat{J}_{y}^{2}+\hat{J}_{z}^{2}$ determine different invariant subspaces. Its maximum value, given by $j=\mathcal{N}/2$, defines the symmetric non-degenerate atomic subspace that includes the ground-state. The Hamiltonian $\hat{H}_{D}$ commutes with the parity operator $\hat{\Pi}=e^{i\pi\hat{\Lambda}}$, where $\hat{\Lambda}=\hat{a}^{\dagger}\hat{a}+\hat{J}_{z}+j\hat{1}$. The operator $\hat{\Lambda}$ represents the total number of excitations with eigenvalues $\lambda=n+m+j$, where $n$ is the number of photons, $m+j$ is the number of excited atoms, and $m$ is the eigenvalue of $\hat{J}_{z}$. In all calculations presented below, we consider the positive parity spectrum of the model.

When the coupling parameter reaches a critical value $\gamma_{c}=\sqrt{\omega\omega_{0}}/2$, a second-order quantum phase transition takes place Hepp and Lieb (1973a, b). The system goes from a normal phase ($\bar{\gamma}<1$ with $\bar{\gamma}=\gamma/\gamma_{c}$), where the ground state has no photons and all the atoms are in their lowest level, to a superradiant phase ($\bar{\gamma}>1$), where the ground state has non-zero expectation values for the number of photons and number of excited atoms.

II.1 Level Statistics and Density of States

The classical limit of the Dicke model can be obtained by using Bloch and Glauber coherent states Bastarrachea-Magnani et al. (2014b, a, 2015); Chávez-Carlos et al. (2016), which allows for the identification of the parameters and energy range that lead to chaos. A main signature of classical chaos in the quantum regime is energy-level repulsion Casati1980 ; Bohigas1984 .

As our case study, we choose $\omega=\omega_{0}$ and $j=100$, and we select a coupling parameter in the superradiant phase, $\bar{\gamma}=2$. For these values, chaos is found at excitation energies above $E\approx-1.6\omega_{0}j$ (see Ref. Chávez-Carlos et al. (2016)). We choose an energy well above this threshold, $E_{c}=-0.5\omega_{0}j$, for which the whole energy shell is covered by chaotic trajectories (see Ref. Chávez-Carlos et al. (2016)). In Fig. 1 (a), we show the level spacing distribution, denoted by $P(s)$, where $s$ is the spacing between nearest-neighboring unfolded energy levels from an energy interval around $E_{c}$. In quantum systems whose classical counterparts are chaotic, the levels are prohibited from crossing and $P(s)$ coincides with the Wigner surmise Guhr et al. (1998), as indeed confirmed in Fig. 1 (a).

With the classical Hamiltonian, it is possible to estimate the energy averaged density of states (DoS), which is given by the expression Bastarrachea-Magnani et al. (2014b),

[TABLE]

where $y_{\pm}=-\bar{\gamma}^{-1}\left(\bar{\gamma}^{-1}\mp\sqrt{2(\epsilon-\epsilon_{gs})}\right)$ and $\epsilon=\frac{E}{\omega_{0}j}$ is the normalized energy. The ground state energy is $\epsilon_{gs}=-1$ for the normal phase, while it is $\epsilon_{gs}=-\dfrac{1}{2}\left(\bar{\gamma}^{2}+\bar{\gamma}^{-2}\right)$ in the superradiant phase. In Fig. 1 (b), we compare the DoS obtained numerically with the expression (2). The agreement is excellent. It is evident from the figure that for $|\epsilon|\leq 1$, the DoS shows a linear dependence on energy, $\nu(E)\propto E$. Our choice of the value of $E_{c}$ for the studies below falls within this region. Notice also that the DoS in Eq. (2) scales linearly with the number of atoms ($j$ appears explicitly in the beginning of the equation), a property that will be useful below to determine the dependence of the timescales of the model on the number of atoms.

III Survival probability and Initial states

The survival probability, $S_{P}(t)$, is a dynamical observable defined as the probability to find an arbitrary initial quantum state $|\Psi(0)\rangle$ at a later time $t$,

[TABLE]

By writing the initial state in terms of the energy eigenbasis, $|\Psi(0)\rangle=\sum_{k}c_{k}|\phi_{k}\rangle$, where $\hat{H}|\phi_{k}\rangle=E_{k}|\phi_{k}\rangle$ and $c_{k}=\langle\phi_{k}|\Psi(0)\rangle$, the survival probability is

[TABLE]

The short-time evolution depends on the energy distribution of the initial state Lerma-Hernández et al. (2018); Torres-Herrera and Santos (2014), while the long-time dynamics is determined by the properties of the spectrum Torres-Herrera and Santos (2017b); Torres-Herrera et al. (2018). The survival probability has been studied in several different contexts, with early works focusing on deviations from exponential behaviors Khalfin (1958); Fonda et al. (1978) and the quantum speed limit Bhattacharyya (1983).

III.1 Initial States

To disentangle the effects of the spectrum from those of the energy components of the initial state in the behavior of the survival probability, we consider ensembles of initial states defined in a given chaotic energy region with components randomly selected, so that

[TABLE]

Above, $r_{k}$ are positive random numbers from an arbitrary probability distribution. For the numerical simulations presented below, we consider an uniform distribution in the interval $[0,1]$ with $n$-th moments $\langle r_{k}^{n}\rangle=1/(n+1)$. The function $f(E)=\rho(E)/\nu(E)$ is used to guarantee that the initial state has a certain selected profile $\rho(E)$, which is achieved by compensating for changes in the density of states. We consider normalized rectangular and Gaussian profiles given respectively by,

[TABLE]

The profiles are centered at the energy $E_{c}$, where we know that chaos dominates the dynamics. The widths of the rectangular and Gaussian profiles are, respectively, $\sigma_{R}$ and $\sigma_{G}$. The lower and upper energy bounds of the Gaussian profile are $E_{\textrm{min}}$ and $E_{\textrm{max}}$, and ${\cal C}$ is a normalization factor,

[TABLE]

with erf being the error function. In the context of quench dynamics, where the system is initially prepared in a coherent state, the energy distribution of the initial state is indeed Gaussian, which makes the Gaussian profile a realistic choice (for some examples, see Ref. Lerma-Hernández et al. (2018)). The bounds $E_{\textrm{min}}$ and $E_{\textrm{max}}$, especially $E_{\textrm{min}}$, are also plausible, since in quantum systems there is always at least a ground state, whose presence should affect the dynamics by partially reconstructing the initial state Muga et al. (2009); Távora et al. (2016, 2017).

In Fig. 2, we show three cases of energy profiles of the initial state, one rectangular and two Gaussian profiles. The numerical results are obtained by averaging over ensembles of 500 initial states. The agreement with the analytical profiles from Eq. (8) and Eq. (11) confirms that 500 is a sufficiently large number to obtain stable results.

III.2 Survival Probability: Before the Correlation Hole

The energy distribution of $|\Psi(0)\rangle$ determines the initial decay of $S_{P}(t)$. This can be seen by expressing the survival probability in Eq. (4) as

[TABLE]

where $\rho_{0}(E)=\sum_{k}|c_{k}|^{2}\delta(E-E_{k})$ is the local density of states (LDoS) or strength function, that is the energy distribution weighted by the components of the initial state. If we approximate the LDoS by its smoothed profile, $\rho_{0}(E)\approx\rho(E)$, we obtain for the rectangular profile,

[TABLE]

and for the Gaussian profile,

[TABLE]

For very short times, $t\ll\sigma_{R,G}^{-1}$, both $S_{P}^{R}(t)$ and $S_{P}^{G}(t)$ show the universal quadratic decay of the survival probability $1-\sigma_{R,G}^{2}t^{2}$. For longer times, both profiles lead to a power-law decay $\propto t^{-2}$. This behavior is evident in Eq. (13) and it can be obtained from Eq. (14) by analyzing it at long times, in which case Távora et al. (2016, 2017),

[TABLE]

where

[TABLE]

Power-law decays of the survival probability are caused by the presence of energy bounds in the LDoS Khalfin (1958) and the power-law exponent depends on how the bounds are approached Urbanowski (2009).

III.3 Survival Probability: Analytical Expression

The expressions in Eqs. (13) and (14) describe accurately the initial decay of $S_{P}(t)$, for which just the shape and bounds of the envelope of $\rho_{0}(E)$ matters. However, the spectra of finite quantum systems are discrete and, in our case, the eigenvalues are correlated. This results in two additional features to the evolution of $S_{P}(t)$, beyond the power-law behavior, which are not captured by Eqs. (13) and (14). They are the manifestations of the spectrum correlations, which appear at long times, and the saturation of the dynamics to the asymptotic value

[TABLE]

around which the survival probability oscillates after relaxation.

To obtain an equation for the full dynamics, we write the survival probability as

[TABLE]

where

[TABLE]

is the so-called inverse participation ratio, which gives the asymptotic temporal value of $S_{P}(t)$. The $I_{PR}$ is a measure of the inverse of the number of elements of a given basis (the energy eigenbasis, in our case) participating in an arbitrary quantum state ($|\Psi(0)\rangle$, in our case).

For the considered ensembles of initial states, it is possible to derive accurate estimates for the ensemble averaged $\langle I_{PR}\rangle$ by using (see Appendix A)

[TABLE]

which, considering random variables $r_{k}$ uniformly distributed, gives for the rectangular profile,

[TABLE]

and for the Gaussian profile,

[TABLE]

Above, $\nu_{c}=\nu(E_{c})$ is the DoS evaluated at the central energy $E_{c}$, which equals the inverse mean spacing of consecutive energy levels in the region probed by the initial state.

Since the components $c_{k}$ of the initial state are random numbers, to compute the ensemble average of $S_{P}(t)$, we can treat the statistical properties of the components and of the spectrum separately. Because the latter has level statistics comparable to that of random matrices from Gaussian orthogonal ensembles (GOE), as shown in Fig. 1 (a), we can follow steps similar to the ones described in Torres-Herrera et al. (2018); Schiulaz et al. to obtain (see Appendix A for details),

[TABLE]

Above, $\eta$ is the effective dimension of the ensemble (see Appendix A) defined as,

[TABLE]

where the second equality is obtained using Eq. (20). In Eq. (23), $S_{P}^{bc}(t)$ describes the behavior of the survival probability before the manifestations of the correlations between the eigenvalues ($bc$ stands for “before correlations”), as given by Eqs. (13) and (14) for the rectangular and Gaussian profiles. This behavior holds until $\left\langle S_{P}(t)\right\rangle$ reaches its minimum value, which is actually below $\left\langle I_{PR}\right\rangle$. Beyond that, the dynamics becomes controlled by the two-level form factor,

[TABLE]

where $\Theta$ is the Heaviside step function. The two-level form factor brings the survival probability from its minimum value up to the asymptotic value, creating the dip that is known as correlation hole Leviandier et al. (1986); Guhr and Weidenmüller (1990); Wilkie and Brumer (1991); Guhr et al. (1998); Alhassid and Levine (1992). The hole is a direct signature of the presence of correlated eigenvalues, and it does not develop in systems with uncorrelated eigenvalues. The equation used above for $b_{2}(\bar{t})$ is the same used for GOE full random matrices Mehta (1991). This implies that beyond the minimum of the correlation hole, the dynamics shows universal properties.

The analytical expression for the survival probability in Eq. (23) describes the complete evolution of $\left\langle S_{P}(t)\right\rangle$, from $t=0$ to saturation. The equation has no fitting parameters. All the parameters entering in Eq. (23) can be determined from the properties of the model and the energy profile of the initial state. As we show in the next section, this analytical expression shows remarkable agreement with our numerical results.

IV Comparing numerical and analytical results

In Fig. 3, we compare numerical results for the survival probability with the analytical expression given by Eq. (23). The light (gray) lines represent the numerical results obtained with a single initial random state for the rectangular (a) and Gaussian [(b) and (c)] energy profiles. The darker (blue) line is obtained by performing ensemble averages over 500 random initial states. The bright (green) curve, following extremely well the ensemble average, is the analytical Eq. (23). As clear from all panels, the ensemble average is needed for the hole to be visible.

The initial decay is determined entirely by the energy profile of the initial states, which is the same for every member of the ensemble. An oscillatory decay modulated by a power law $\propto t^{-\lambda}$ is seen in Figs. 3 (a) and (b). As mentioned before, this behavior is caused by the bounds in the energy profiles, which determine also the size of the oscillations. For the rectangular profile (a), the exponent is indeed $\lambda=2$, as obtained for $S_{P}^{bc}(t)$ in Eq. (13), whereas for the strongly bounded Gaussian profile (b), the power-law exponent obtained numerically is $\lambda=1.7$, instead of 2 as in Eq. (15). This is because the oscillations in Fig. 3 (b) start at a temporal scale where the two-level form factor $b_{2}$ is not negligible, so it affects the exponent. The power-law decay in Figs. 3 (a) and (b) is followed by the correlation hole. In the case of the weakly bounded Gaussian profile of panel (c), no trace of the modulated oscillations is left, because they occur at a temporal scale when $S_{P}^{bc}(t)$ is already extremely small with respect to the $b_{2}$ term. In all panels, once saturation is reached, the survival probability only fluctuates around its asymptotic value.

The minimum of the correlation hole is indicated in Figs. 3 (a), (b) and (c) with the lowest horizontal dashed line. The dynamics beyond this minimum point depends on level statistics, as confirmed by the fact that the behavior of the ensemble averages is very well described by the analytical expression (23), where the same two-level form factor $b_{2}$ used for full random matrices was employed. However, the size of the temporal fluctuations after the mimimum depends on the fine details of the particular spectrum and on the level of delocalization of the initial state written in the energy eigenbasis Reimann (2008); Short (2011); Short and Farrelly (2012); Zangara et al. (2013).

Contrary to ensembles of random matrices or disordered models, where one can further reduce the temporal fluctuations of $\left\langle S_{P}(t)\right\rangle$ with averages over many energy spectrum realizations, in the case of the Dicke model, the spectrum is exactly the same for every member of the ensemble. Thus, to further reduce the fluctuations in the ensemble averaged $\left\langle S_{P}(t)\right\rangle$, we perform an additional time average over temporal windows of constant size in logarithmic scale, i.e. temporal windows whose sizes increase exponentially in time. By plotting this temporal averages against the mean time of the respective windows, we obtain the results shown in Figs. 3 (d). This smoothing procedure results in numerical curves that are almost identical to the analytical curves, further validating Eq. (23) and the approach that led to it.

The fact that the dynamics beyond the minimum of the correlation hole is governed entirely by the two-level form factor implies that the time to reach saturation depends only on how the $b_{2}$ function approaches $\left\langle I_{PR}\right\rangle$ (indicated in Fig. 3 with the highest horizontal dashed line). Provided the initial state is fully extended in the energy eigenbasis, counting with the participation of all (most) energy levels in the energy interval characterizing $|\Psi(0)\rangle$, the relaxation time is independent of the initial state. Indeed, as seen in Fig. 3 (d), the time to reach $\left\langle I_{PR}\right\rangle$, which is shown with the rightmost vertical dashed line, is the same for the three different energy profiles.

V Timescales of the survival probability

In hands of the analytical expression for the survival probability, we can derive analytically the timescales involved in the relaxation process. We focus on the two longest timescales: the time to reach the minimum of the correlation hole, referred to as Thouless time $t_{Th}$, and the final relaxation time $t_{r}$.

V.1 Thouless Time

The Thouless time divides the dynamics of chaotic systems in two temporal regions, before $t_{Th}$ the dynamics is governed by the shape of the energy distribution of the initial state and the energy bounds, whereas after $t_{Th}$ the dynamics becomes comparable to that obtained with ensembles of full random matrices. The Thouless time marks the point where the term $\eta S_{P}^{bc}(t)$ in Eq. (23) meets the function $b_{2}(t/2\pi\nu_{c})$, being derived from

[TABLE]

Therefore, we need to examine $S_{P}^{bc}(t)$ at long times and $b_{2}(t/2\pi\nu_{c})$ at short times, i.e. in the temporal range $\sigma^{-1}\ll t\ll\nu_{c}$ [recall that $\sigma$ is the width of the LDoS and $\nu_{c}$ is the inverse of the mean level spacing for the eigenvalues involved in the evolution of $|\Psi(0)\rangle$].

At long times, $\eta S_{P}^{bc}(t)$ shows a power-law decay $\propto t^{-2}$. More specifically, for the rectangular energy profile, the temporal average of the oscillatory decay in Eq. (13) leads to

[TABLE]

For the Gaussian profile, associated with Eq. (14), if the conditions

[TABLE]

are fulfilled, the form of $S_{P}^{bc}(t)$ in the timescale where it meets $b_{2}$ is given by [see Eq. (15)]

[TABLE]

This is what happens for the strongly bounded Gaussian energy profile of Fig. 2 (b). Otherwise, if the conditions in Eq. (28) are not fulfilled, the early Gaussian decay still persists at the meeting point with the $b_{2}$ function and

[TABLE]

This is what happens for the weakly bounded Gaussian energy profile of Fig. 2 (c).

At short times, the two-level form factor is dominated by a linear term,

[TABLE]

Combining Eqs. (27) and (31), we obtain the Thouless time for the rectangular ensemble

[TABLE]

where we used that $\nu_{c}=\frac{\eta}{2\sigma_{R}}$ (see Eq. (24) and Appendix A). From Eq. (29) and Eq. (31), we arrive at the Thouless time for the strongly bounded Gaussian energy profile,

[TABLE]

and if the conditions (28) are not fulfilled, using Eq. (30), we have

[TABLE]

The Thouless time obtained in the equations above are indicated in Fig. 3, showing excellent agreement with the numerics.

V.2 Relaxation Time

The relaxation time depends only on the $b_{2}$ function at long times, which grows toward saturation following a power-law behavior,

[TABLE]

Even though Eq. (27) and Eq. (29) decay with the same power-law exponent 2, as in Eq. (35), the latter is proportional to $\eta^{2}$, since $\nu_{c}\propto\eta$, while the former equations are proportional to $\eta$, which justifies considering only Eq. (35). We define the relaxation time according to

[TABLE]

where $\delta$ is a small parameter determining the point where $\left\langle S_{P}(t)\right\rangle$ is already within the fluctuations around the asymptotic value. We arrive at

[TABLE]

which holds for the three energy profiles. This time is proportional to the inverse of the mean level spacing, $\nu_{c}$, which is the largest timescale of a quantum system and is known as the Heisenberg time. In Fig. 3, $t_{r}$ is indicated with the rightmost vertical lines, showing a excellent agreement with the numerical results.

V.3 Scaling with system size: Thouless and relaxation time

With Eqs. (32), (33), (34) and (37), we can determine the dependence of the Thouless and relaxation times on the size of the system, i.e on the number $\mathcal{N}=2j$ of two-level atoms. Both times depend on $\nu_{c}$, which scales linearly with $j$. We can write $\nu_{c}=\nu_{o}j/\omega$ and evaluate $\nu_{o}$ numerically, which for the chosen energy $E_{c}$ is $\nu_{o}=0.6027$.

The Thouless time depends additionally on the widths $\sigma_{R}$ and $\sigma_{G}$ of the energy distribution of the initial state and, for the Gaussian profile, on the energy bounds $E_{min}$ and $E_{max}$. The scaling $\sigma_{R,G}=\sigma_{R,G}^{o}j^{\beta}$ of these widths, as well as the scalings of $(E_{c}-E_{min})\propto j^{\alpha_{1}}$ and $(E_{max}-E_{c})\propto j^{\alpha_{2}}$, can in principle be selected at will in the range $-1\leq\beta\leq 1$ and $-1\leq\alpha_{i}\leq 1$. The lower values $\beta=\alpha_{i}=-1$ are imposed by the scaling of the mean-level spacing of consecutive energy levels and the upper bound is given by the scaling of the energy spectrum, which is proportional to $j$. A physical relevant choice for the previous scalings is $\beta=1/2$, which is the scaling of the energy widths of minimal uncertainty coherent states Schliemann (2015), and $\alpha_{i}=1$, which implies that the bounding energy interval of the Gaussian profile scales as the energy spectrum. Therefore for the rectangular profile and strongly bounded Gaussian profile satisfying conditions (28), the Thouless time scales as

[TABLE]

where $t_{o}^{R,G}$ is a constant determined by $\nu_{o}$, and $\sigma_{R,G}^{o}$. For the rectangular profile, this scaling is valid for any $j$, but for the Gaussian profile, it is valid up to a finite value of $j$. This is because we assume that $(E_{c}-E_{min})$ and $(E_{max}-E_{c})$ grow with $j$ faster than $\sigma_{G}$, which implies that for large enough $j$ the conditions (28) will not be satisfied anymore, switching to the scenario of Eq. (34). For the weakly bounded Gaussian profile, described by Eq. (34), the Thouless time for large $j$ is given by,

[TABLE]

where $c_{o}=4\pi^{(3/2)}(\sigma_{G}^{o})^{2}\nu_{o}^{2}/\omega^{2}$, and we have approximated the error functions by their asymptotic values, $\lim_{z\rightarrow\infty}\textrm{erf}(z)\rightarrow 1$.

The relaxation time is independent of the details of the initial state and scales linearly with $\nu_{c}$, so it is given simply by

[TABLE]

The distance between the Thouless and the relaxation time diverge with $j$, which means that the correlation hole gets elongated as the number of atoms increases.

V.4 Depth of the correlation hole

With the Thouless time, we can quantify the relative depth of the correlation hole through the expression

[TABLE]

For the ensembles considered in this paper, we can calculate the depth of the correlation hole for $j\gg 1$. A direct substitution of the Thouless time in the analytical expression for the survival probability in Eq. (23) allows to demonstrate that

[TABLE]

where in the last step we have used Eq. (24). From this result, we obtain

[TABLE]

The value $\kappa_{\infty}$ gives an upper bound for the depth of the correlation hole for finite $j$. In the case of random variables $r_{k}$ uniformly distributed in the interval $[0,1]$, as considered here noteGOE , this bound is $\kappa_{\infty}=3/4$.

The actual values of $\kappa$ for the finite systems studied, where $j=100$, are obtained by substituting Eqs. (32), (33) and (34) in Eq. (23), which gives $\left\langle S_{P}(t_{Th})\right\rangle$, and by getting $\left\langle I_{PR}\right\rangle$ from Eqs. (21) and (22). We get $\kappa=0.672$ for the rectangular ensemble, $\kappa=0.711$ for the ensemble from the strongly bounded Gaussian profile, and $\kappa=0.748$ for the ensemble from the weakly bounded Gaussian profile. The analytical estimates for $\left\langle S_{P}(t_{Th})\right\rangle$ and $\left\langle I_{PR}\right\rangle$ are depicted in Fig. 3 (a), (b), and (c) with horizontal lines, showing excellent agreement with the numerical results. This confirms that the analytical expression in Eq. (23) describes the survival probability at any timescale.

VI CONCLUSIONS

We obtained an analytical expression that describes remarkably well the entire evolution of the averaged survival probability, $\langle S_{P}(t)\rangle$, for the Dicke model in the chaotic regime and allowed us to derive analytical expressions for the different timescales involved in the relaxation to equilibrium. Due to spectral correlations, the survival probability exhibits a correlation hole. We find that the initial decay of $\langle S_{P}(t)\rangle$ and the time $t_{Th}$ for it to reach the minimum of the correlation hole (Thouless time) depend on the energy profile of the initial states. Beyond the Thouless time, the dynamics is universal, being governed by the two-level form factor of the GOE. This implies that the time, $t_{r}$, for the survival probability to reach equilibrium (relaxation time) depends only on the energies, being proportional to the inverse of the mean level spacing. An interesting future direction would be to investigate the timescales involved in the relaxation process of other physical observables that are relevant for the Dicke model.

Acknowledgements

We acknowledge the support of the Computation Center- ICN, in particular to Enrique Palacios, Luciano Diaz and Eduardo Murrieta. D.V, J.G.H, M.A.B.-M. and S.L.-H acknowledge Jorge Chávez-Carlos for fruitful discussions and his technical support. E.J.T.-H. acknowledges funding from VIEP-BUAP (Grant Nos. MEBJ-EXC19-G, LUAG- EXC19-G), Mexico. L.F.S. is supported by the NSF Grant No. DMR-1603418. S.L.-H. acknowledges financial support from Mexican CONACyT project CB2015-01/255702, J.G.H. and D.V. acknowledge funding from DGAPA- UNAM project IN109417.

Appendix A Ensemble averages of the survival probability

In this appendix we present the steps involved in the derivation of the analytical expression given by Eq. (23) for the ensemble average of the survival probability. We begin with Eq. (18) and perform ensemble averages, taking into account that the eigenvalues and the components $|c_{k}|^{2}$ of the initial state are statistically independent,

[TABLE]

Let us consider first the ensemble average of $I_{PR}$. Using Eq. (5), we have

[TABLE]

For a large number of components, the average of the second line above can be approximated as

[TABLE]

where in the last equality we have used the fact that $\langle r_{q}^{n}\rangle$ is actually independent of index $q$. By inserting this result in Eq. (44), we obtain

[TABLE]

Here, we have introduced the effective dimension of the ensemble

[TABLE]

whose name comes from the fact that it reduces to the number of states participating in the rectangular ensemble, as it is shown below. We now approximate the sums in Eq. (46) by integrals,

[TABLE]

to obtain

[TABLE]

where we have used $f(E)=\rho(E)/\nu(E)$, and, in the last equality, the normalization of $\rho(E)$. Finally, since $\nu(E)$ varies linearly in the energy interval where $\rho(E)$ is significant, we substitute the function $\nu(E)$ by its value in the center of the profile distribution $\nu_{c}\equiv\nu(E_{c})$ and obtain

[TABLE]

and

[TABLE]

From the expression for $\eta$, it is clear that, in the case of the rectangular profile $\eta=2\nu_{c}\sigma_{R}$, which gives the number of states participating in the ensemble.

For the first term in Eq (43), we have to evaluate the ensemble average

[TABLE]

where $(l\not=k)$ and, for simplicity, we introduced the shorthand notation $f_{k}\equiv f(E_{k})$. To obtain an approximation to the average, we consider the identity,

[TABLE]

From this, we obtain

[TABLE]

By taking the ensemble average of this expression and assuming that

[TABLE]

where $A$ is a constant independent of indexes $l$ and $k$, we get

[TABLE]

which implies that

[TABLE]

The sum in the denominator of the equation above can be expressed in terms of the effective dimension $\eta$ [see Eq. (47)], as follows

[TABLE]

where in the last step, the normalization $\sum_{q}f_{q}\rightarrow\int dE\rho(E)=1$ was used. With the the above result, the average can be written as

[TABLE]

which, when substituted in Eq. (43), leads to

[TABLE]

We now turn our attention to the double sum $\sum_{k\neq l}$ in Eq. (49). To solve it, we use

[TABLE]

where the Dyson two-point correlation function, $R(E,E^{\prime})=\nu(E)\nu(E^{\prime})-T(E-E^{\prime})$, includes the DoS, $\nu(E)$, and the two-level cluster function, $T(E-E^{\prime})$, which has information about the correlations between the eigenvalues Mehta (1991). We then obtain

[TABLE]

Using unfolded energy variables, which leads to universal functions in the limit of an infinite number of levels, we have Mehta (1991)

[TABLE]

With this function, the second integral in Eq. (50) is

[TABLE]

which, in terms of variables $z=E^{\prime}$ and $x=(E-E^{\prime})\nu_{c}$ becomes

[TABLE]

with $\tilde{t}=t/(2\pi\nu_{c})$. By expanding $\rho(z+x/\nu_{c})$ in powers of $x$ and considering only the lowest order, the double integral can be approximated by a product of two independent integrals

[TABLE]

where we have used the effective dimension introduced before, and $b_{2}(\tilde{t})$ is the known Fourier transform of the GOE-two level cluster function Mehta (1991), the so called two-level form factor shown in Eq. (25). We use the same $b_{2}(\tilde{t})$ as in GOE matrices, because the unfolded spectrum of the Dicke model has correlations comparable to those of the GOE levels.

Gathering the previous results together, we obtain for the ensemble average of the survival probability

[TABLE]

Since the squared absolute value inside the parenthesis is the short time decay $S_{p}^{bc}$ given by Eqs. (13) and (14), we finally arrive at our analytical expression in Eq. (23).