Unifying scrambling, thermalization and entanglement through measurement of fidelity out-of-time-order correlators in the Dicke model.

Scrambling is the process by which information stored in local degrees of freedom spreads over the many-body degrees of freedom of a quantum system, becoming inaccessible to local probes and apparently lost. Scrambling and entanglement can reconcile seemingly unrelated behaviors including thermalization of isolated quantum systems and information loss in black holes. Here, we demonstrate that fidelity out-of-time-order correlators (FOTOCs) can elucidate connections between scrambling, entanglement, ergodicity and quantum chaos (butterfly effect). We compute FOTOCs for the paradigmatic Dicke model, and show they can measure subsystem Rényi entropies and inform about quantum thermalization. Moreover, we illustrate why FOTOCs give access to a simple relation between quantum and classical Lyapunov exponents in a chaotic system without finite-size effects. Our results open a path to experimental use FOTOCs to explore scrambling, bounds on quantum information processing and investigation of black hole analogs in controllable quantum systems.

Introduction

Recent studies have shown that isolated many-body quantum systems, under unitary time evolution, can become highly entangled and thus thermalize. This understanding has led to insights as to how statistical mechanics emerges in closed quantum systems[1-3]. Moreover, the relevance of entanglement as a resource for quantum information processing, quantum communication and metrology has stimulated cross-disciplinary efforts to quantify and characterize entanglement. Experimental progress in controlling clean, highly isolated, and fully tunable quantum systems, where entanglement can be measured, have resulted in radical advances in this direction. However, such measurements have been restricted to few body systems, including arrays of 6 × 2 bosonic atoms[4], three superconducting qubits[5], and systems of trapped ions[6,7]. The model we study here and the measurements we propose can be implemented in trapped ions with more than 100 spins.

Concurrently, out-of-time-order correlations (OTOCs)[8-13]have been identified as measures of the dynamics of quantum information scrambling. Here, , with a quantum many-body Hamiltonian, and and two initially commuting and unitary operators. While OTOCs can be computed with respect to any (possibly mixed) state, here we focus on the case where the initial state of the system is pure. The quantity encapsulates the degree that and fail to commute at later times due to the time evolution of under . (If is not unitary but a projector, e.g. and commutes with the density matrix of the initial state, then .) The fastest scramblers[8-10,14], such as black holes, feature an exponential growth of scrambling which manifests as . Here, λQ is the quantum Lyapunov exponent that serves as a proxy for quantum chaos. Regardless of the OTOCs’ apparent complexity[13,15-17], the capability to perform many-body echoes (see Fig. 1) in current experiments[18-21] has opened a path for the experimental investigation of quantum scrambling; however, so far those have not probed quantum chaos or fast scrambling.

Fig. 1

Unifying chaos, scrambling, entanglement and thermalization through the measurement of fidelity out-of-time-order correlators (FOTOCs). a Scheme: an initial state, |ψ0〉 is evolved under an interacting Hamiltonian for a time t. Inverting the sign of and evolving again for time t to the final state |ψf〉, implements the many-body time-reversal, which ideally takes the system back to the initial state |ψ0〉. If a perturbation is inserted between the two halves of the time evolution and the many-body overlap with the initial state is measured at the end of the protocol, , then a special type of fidelity OTOC (FOTOC) is implemented. b The Dicke model is engineered in a Penning trap ion crystal by applying a pair of lasers, resonant only with the center-of-mass mode, to generate the spin−phonon interaction and resonant microwaves to generate the transverse field

Here we show that fidelity out-of-time-order correlators (FOTOCs), a specific family of fidelity out-of-time-order correlators, which set to be a projector on the initial state, can provide profound insight on scrambling behavior. We explicitly compute FOTOCs in the in the Dicke model[22], an iconic model in quantum optics, and illustrate how FOTOCs elucidate theoretical connections between scrambling, volume-law Rényi entropy (RE) and thermalization, while linking quantum and classical chaos (Fig. 1a). Additionally, we discuss how one can probe these connections readily in experiments.

Results

Model

The Dicke model (DM)[22] describes the coupling of a single large spin and a harmonic oscillator and has been recently implemented in atomic[23-26] and trapped ion setups[27]. The Hamiltonian of the DM is given bywhere B characterizes the strength of the transverse field, δ the detuning of the bosonic mode from the driving field with strength g that generates the spin−boson coupling. Here, g, δ, B ≥ 0. The operator () is the bosonic annihilation (creation) operator of the mode, and are collective spin operators with (α = x, y, z) the Pauli matrices for the jth spin-1/2.

Connections between scrambling dynamics and chaos

Even when restricted to the Dicke manifold, i.e. states with S = N/2, with S(S + 1) the eigenvalue of the total spin operator , this model exhibits rich physics (see Fig. 2a). At zero temperature, T = 0, the DM features a quantum phase transition (QPT) as the system crosses a critical field Bc = 4g2/δ. For B > Bc (normal phase), the ground-state is described by spins aligned along the transverse field and a bosonic vacuum. For B < Bc (superradiant phase), the ground-state is ferromagnetic, , and characterized by macroscopic occupation of the bosonic mode (Fig. 2a). Furthermore, in the superradiant phase (B < Bc), the DM features a family of excited-state quantum phase transitions (ESQPTs). The ESQPTs are signaled by singularities in the energy-level structure and a change in the spectral statistics[28-31] at a critical energy Ec = −BN/2 that coincides with the ground-state energy of the normal phase. Figure 2a shows how the nearest-neighbor spacing distribution P(s), where s is a normalized distance between two neighboring energy levels, features a different character on either side of Ec. For E > Ec the spectral statistics are similar to the Wigner−Dyson distribution , which in random-matrix theory describes a chaotic system. For E < Ec the shape of the histograms is neither Wigner−Dyson nor Poissonian PP(s) = exp(−s). The latter characterizes level statistics of non-ergodic systems, and is observed in the normal phase. While the deviations from clear Wigner−Dyson or Poissonian statistics in regimes II and III are attributable to finite-size effects[28], we emphasize that even for this small system they clearly show a stark contrast in the degree of level repulsion, which is a qualitative signature of quantum chaos.

Fig. 2

Characterization of classical and quantum chaos in the Dicke model. a Phase diagram of the Dicke model. At zero temperature it exhibits a quantum phase transition between a normal to a superradiant phase, at B = Bc. A line of excited energy quantum phase transitions (ESQPTs) occurs at the critical energy Ec = −BN/2, signaled by singularities in the energy level structure (indicated by stars). Note that for figure clarity we have used a small system N = 20 resulting in the small deviation of the ESQPTs from Ec = −BN/2. The ESQPTs are accompanied by a change in the level statistics which we denote by (I)−(III) (note that no eigenstates exist in the unlabeled white region). For (II) and (III) the spectrum is divided into low and high energy parts, separated by the ESQPT at E = Ec, from which the statistics P(s), where s is the level spacing, are computed separately. (I) exhibits Poissonian statistics (regular regime), while (II) displays statistics similar to a Wigner−Dyson distribution indicative of level repulsion and quantum chaos, and (III) exhibits a mixture of both. The numerical parameters are g/(2π) = 0.66 kHz and δ/(2π) = 0.5 kHz. b Lyapunov exponents for the mean-field dynamics of an ensemble of random states sorted by normalized mean-field energy E/|Ec| with Ec = −BN/2, as a function of the field relative to critical field Bc = 4g2/δ. A crossover between regular (B > Bc) and chaotic dynamics (B < Bc) characterized by and λL > 0 respectively, occurs at B = Bc. For B > Bc and energies the dynamics becomes increasingly regular. Source data are provided as a Source Data file

Similar features appear in the classical dynamics of the DM[30,32-35], manifested in the different behavior of trajectories in phase-space computed from the mean-field equations of motion for: , where denotes the expectation values, and αR(I) is the real (imaginary) part of . In the superradiant phase and for mean-field energies , two trajectories initially separated by in phase-space diverge as at sufficiently long times[36]. The exponential growth, associated with a positive Lyapunov exponent λL > 0, diagnoses chaos in a classical system. In Fig. 2b we show the maximal Lyapunov exponent for an ensemble of random initial product states as a function of the transverse field and the normalized mean-field energy E/Ec (see Methods). For E < Ec in the superradiant phase (B < Bc) and all energies in the normal phase (B > Bc), the Lyapunov exponent is small or zero, consistent with the Poissonian character of the quantum-level statistics in this parameter regime[34,35]. For E > Ec and B < Bc a positive exponent is found signaling chaos. Note that the state , where , lies exactly at the ESQPT critical energy, , and possesses the largest classical λL (see Fig. 2).

In quantum systems OTOCs may serve as a diagnostic for quantum chaos. However, such diagnosis has proved difficult, since any exact numerical treatment is only possible in small systems, where many-body observables saturate quickly at the Ehrenfest time given by , at which the quantum information is thoroughly lost to a “local” observer. Here we demonstrate that we can overcome this limitation and compute OTOCs for macroscopic systems if, for a Hermitian operator , one restricts to be a sufficiently small perturbation and sets to be a projection operator onto a simple initial state |ψ0〉, i.e. . This is because in the perturbative limit , this particular type of fidelity OTOC (FOTOC)[18,19], (such that for a pure state ) reduces to[37]

where is the variance of . This relation establishes a connection between the exponential growth of quantum variances and quantum chaos, enables us to visualize the scrambling dynamics of a quantum system using a semi-classical picture[38] and to map the FOTOC to a two-point correlator which can be computed using well-known phase-space methods, such as the truncated Wigner approximation (see Methods)[39,40]. We observe perfect agreement between the exact dynamics of the FOTOC with the associated variance, for sufficient small δϕ, enabling us to use phase-space methods to compute the FOTOCs in a parameter regime inaccessible to exact numerical diagonalization where exponential scrambling can be clearly identified.

Moreover, it provides a link between the FOTOCs and the quantum Fisher information (QFI)[19,41-43], as the variance of is proportional to the QFI of a pure state, whilst for a mixed state the variance gives a lower bound on the QFI. Note that in the latter case FOTOCs are defined by replacing by the initial density matrix , and expectation values are computed by appropriate traces. The QFI quantifies the maximal precision with which a parameter δϕ in the unitary of can be estimated using an interferometric protocol with an input quantum state |ψ(t)〉, while simultaneously serving as a witness to multipartite entanglement[19,44-46].

In Fig. 3a we plot the FOTOCs of a small perturbation using starting with . In the superradiant phase we observe that after a short time of slow dynamics, , the FOTOCs feature an exponential growth , before saturating at (see inset). The quantum exponent is found to be independent of system size N. For this initial state, and all the product states we have investigated numerically (Supplementary Note 1), we have observed that , as shown in Fig. 3b. Indeed, for any that corresponds to a linear function of the classical phase-space variables (see Methods and Supplementary Note 1), the quantum exponent should be related to the classical Lyapunov exponent by this relation. A similar factor of two relating the classical and quantum exponents has previously been observed in refs. [37,47-49]. This correspondence can be explained by semi-classical arguments (see Methods), and the numeric prefactor is attributable to the definition of the classical Lyapunov exponent in terms of a distance in phase-space, while the FOTOC reduces to the quantum variance.

Fig. 3

Signatures of classical chaos in quantum FOTOCs. a Initial exponential growth of the FOTOC, and the initial state (see Supplementary Note 1 for examples of exponential growth in other states). We assume such that we may equivalently use for the plotted data. The scrambling time t* is defined by the saturation of the FOTOC, which we extract from the first maximum and plot in the inset (blue data). We find with a0 a fit parameter (gray line). b Lyapunov exponent, λ, as a function of transverse field: Quantum λQ (red markers) and classical 2λL (solid lines). Superscript notation of the exponents denotes the initial polarization of the chosen coherent spin state. Top panel for , the same state as (a), and bottom for , here N = 104 particles. In both plots we observe . Error bars for λQ are a 95% confidence interval from an exponential fitted to the numerical data. Coupling g and detuning δ are same as Fig. 2. In (a) B/(2π) = 0.7 kHz (B/Bc = 0.2). Source data are provided as a Source Data file

FOTOCs as a probe of entanglement and quantum thermalization

We now move beyond the semi-classical arena and explore connections between FOTOCs and entanglement entropy. In a closed system S the second-order Rényi entropy, RE, measures the entanglement between a subsystem A and its complement Ac = S − A, with the reduced density matrix of A after tracing over Ac. Although scrambling and entanglement buildup are closely connected, they are not the same. Nevertheless, a formal relationship between the OTOCs and exists[11], which requires averaging OTOCs over a complete basis of operators of the system subsystem A. Based on this relation, measuring RE via OTOCs appears as challenging as directly measuring . However, this is not always the case. We will show that for collective Hamiltonians, such as the DM, there is a simple correspondence between the Fourier spectrum of FOTOCs and the RE, which facilitates experimental access to via global measurements and collective rotations.

To illustrate the connection we first write the density matrix of the full system in a basis spanned by the eigenstates of the spin operator where e is a unit vector in the Bloch sphere, satisfying , and the mode number operator , i.e. . We adopt a convention for the coefficients of the density matrix elements where superscripts are associated with the bosonic mode, and subscripts with the spin. In this basis the density matrix can be divided into blocks, with , in such a way that contains all coherences between states with spin eigenvalues that differ by M. A similar decomposition can be performed in terms of the bosonic coherences as with . Associated with this representation one can define the so-called multiple quantum intensities . Of particular interest for us are the components which are “incoherent” with respect to .

The intensities can be accessed experimentally from FOTOCs via the relation [18,19,50-52] by choosing and or , i.e. collective spin or boson rotations respectively. In terms of the the entanglement between the spins and the phonons characterized by the purity can be written as

The terms and are explicitly detailed in the Methods, but importantly is composed purely of the diagonal elements of while contains information about coherences. During unitary evolution the characteristic dephasing time of the coherences is , which for scrambling systems is much faster than . After tc any remaining coherences are fully randomized and destructively interfere yielding . This feature, together with the fact that for those systems also the magnitude of becomes much smaller than and as the density matrix spreads out over the systems degrees of freedom, allows us to approximate . While at t < tc these conditions are not necessarily satisfied, we still find that there can be a correspondence between the FOTOCs and RE by picking a state that is fully incoherent at time t = 0, . An example of such a state is and . This choice enforces the term to remain small at short times. Moreover, for |Ψ0〉 we find it is also possible to access via even in the regime B > Bc, where no scrambling occurs. This is because the contributions from and cancel and .

In Fig. 4a we show the typical behavior of the RE, , in the two different phases for |Ψ0〉. First, in the normal phase (panel (i)), B ≪ Bc, the dynamics is dominated by precession about the transverse field and the entanglement entropy exhibits small amplitude oscillations[53]. Conversely, in the superradiant phase (panel (ii)) B ≪ Bc we observe a rapid growth of entanglement and saturation past the transient regime. We summarize our results in Fig. 4b where we plot the time-averaged value of vs. . We associate the fast growth of at with a crossover from the integrable to the chaotic regime. To further illustrate this connection, we compare the approximate RE obtained via and with the exact RE in Fig. 4b. It is observed that in all parameter regimes one can make a quantitative link between the RE and FOTOCs, especially under proper optimization of the rotation axis at each time to minimize the coherence and diagonal terms in Eq. (4) (see Methods and Supplementary Methods).

Fig. 4

Using RE and FOTOCs to characterize chaos and thermalization in the Dicke model. a Time evolution of the spin−phonon RE (black lines) for the initial state with B > Bc (top) and B < Bc (bottom). The RE is tracked excellently by the FOTOC expressions (blue lines) and respectively. Here, is chosen to minimize the coherence and diagonal terms in Eq. (4) (Supplementary Methods). b Long-time spin−phonon RE as a function of transverse field. To remove finite-size effects and residual oscillations we plot a time-averaged value for 4 ms ≤  t ≤ 12 ms (FOTOC quantities are averaged identically). The regular and chaotic dynamics for the initial state are clearly delineated: for B > Bc and for B < Bc respectively. Error bars indicate standard deviation of temporal fluctuations. In the inset we plot the same FOTOC quantities but including decoherence due to single-particle dephasing at the rate Γ = 60 s−1. The coherent parameters g, B and δ are enhanced by a factor of 16 compared to the main panel, as per ref. [60]. c Time-averaged distribution functions (markers) for spin-projection P(Mz) and phonon occupation P(n) (6 ms ≤ t ≤ 12 ms). We compare to the distribution of the diagonal ensemble (purple bars, see Methods). d Bipartite RE (black markers) as a function of partition size L of the spins, averaged over same time window as (c). For comparison, we plot the RE of a thermal canonical ensemble with corresponding temperature T fixed by the energy of the initial state , and the RE of the diagonal ensemble (see Methods). Volume-law behavior of the RE is replicated by the FOTOC quantity (blue markers). Note that the dimension of the spin Hilbert space scales linearly with L. Shaded regions indicate standard deviation of temporal fluctuations. Data for (a)–(d) is obtained for N = 40, with g and δ identical to calculations of Fig. 2. For (c) and (d) we choose B/(2π) = 0.7 kHz (B/Bc = 0.2). Source data are provided as a Source Data file

The saturation of for B < Bc is a signature of thermalization. One can test how “thermalized” the quantum system is by comparing the behavior of the spin and phonon distributions in the long time limit with those of the corresponding diagonal ensemble, characterized by a mixed density matrix with purely diagonal elements (see Methods)[1-3]. These comparisons are shown in Fig. 4c, where the time evolved distributions and the ones drawn from the diagonal ensemble are almost indistinguishable.

We can also investigate the growth of entanglement on different size bipartitions for B < Bc. For that we split the spin system into a subsystem of size L ≤ N and evaluate by computing the reduced density matrix by tracing over the bosonic degree of freedom and the remaining N − L spins. To demonstrate the entanglement grows with system size in a manner consistent with an equivalent thermal state we plot the predictions of a canonical ensemble (see Methods). We observe volume-law entanglement growth for L ≪ N (see Fig. 4d). However, for L ~ N the entanglement growth deviates from this simple prediction. These deviations occur as the full state of the system is pure, and thus eventually one needs to recover , requiring a negative curvature. To demonstrate the intertwined nature of thermalization and the buildup of entanglement we plot the predictions of a canonical ensemble indicated by the dotted purple line (see Methods). We note that FOTOCs can also be used to probe this scaling of the RE with subsystem size. To this end, both and should be restricted to a partition of size L of the system, but otherwise the corresponding multiple quantum intensities are computed as discussed above (see also Methods). Figure 4d shows the excellent agreement between the partial system FOTOCs (blue squares) and RE (black diamonds), comparisons that illustrate the utility of FOTOCs to characterize complex many-body entanglement.

Experimental implementation in trapped ion simulators

Trapped ions present a promising experimental platform for the investigation of the physics discussed here[27,54,55]. Here we focus on two-dimensional arrays in a Penning trap where a tunable coupling between the ion’s spin, encoded in two hyperfine states, and the phononic center-of-mass (COM) mode of the crystal can be implemented by a pair of lasers with a beatnote frequency detuned by δ from the COM mode and far from resonance to all other modes, which remain unexcited (Fig. 1b). In the presence of microwaves (which generate the transverse field) resonant with the spin level splitting, the effective Hamiltonian is of the form of Eq. 2 as benchmarked in refs. [27,56]. The dynamical control of the transverse field and sign of the detuning from the COM mode enables straightforward implementation of a time-reversal protocol to measure FOTOCs[19] (see Fig. 1). Additionally, the many-body echo requires the application of a spin echo π pulse along which reverses the signs of and simultaneously.

Our proposal requires the ability for measuring the fidelity of the full spin−phonon state, which we have not yet demonstrated experimentally. However, this will be possible through a generalization of the protocol discussed in ref. [57] (see Methods). Additionally, our proposal can be adversely affected by decoherence present in the experiment. However, the impact of decoherence will be minimized in future experiments by increasing the magnitude of relevant couplings of the DM via parametric amplification of the ions’ motion[27,58], thus reducing the ratio of dissipative to coherent evolution. We illustrate the predicted effect of decoherence, which is dominated by single-particle dephasing due to light scattering from the lasers, in the inset of Fig. 4b. We include the enhancement of the coherent parameters via the protocol described in ref. [58] while using the typical experimental decoherence rate of Γ = 60 s−1. The single-particle decoherence is modeled by an exponential decay of the FOTOC components (see Methods). The numerical calculation indicates that even with decoherence the crossover between the two regimes at B ~ Bc is still well captured. Due to numerical complexity of solving a master equation we restrict our simulations to N = 40 ions.

Discussion

We have demonstrated that FOTOCs connect the fundamental concepts of scrambling, chaos, quantum thermalization, and multipartite entanglement in the DM. While the concepts presented here have been limited to collective Hamiltonians, we believe they can be generalized to more complex many-body models (Supplementary Note 2). For example, FOTOCs could provide an alternative approach for performing efficient measurements of RE in a way comparable to other state-of-the-art methods which have been used to probe entanglement in systems with up to 20 ions[7]. Generically, FOTOCs could serve as an experimental tool capable of uncovering bounds on information transport and computational complexity, and shed light on how classical behaviors in macroscopic systems emerge from purely microscopic quantum effects.

Note added: Upon completion of this manuscript we became aware of the recent preprints[59,60], which present the numerical and analytic investigation of OTOCs in the Dicke model.

Methods

Classical dynamics and equations of motion

The results presented for the classical model in Fig. 3 are obtained from the Heisenberg equations of motion for the operators via a mean-field ansatz, wherein the operators are replaced by the c-number expectation values, i.e., for and (where we adopt as the real (imaginary) component of ). We thus obtain an equation of motion for ,

where

Lyapunov exponent

The existence of classical chaos can be characterized by the Lyapunov exponent λL. By definition, classical chaos implies that two initially close trajectories separated by a distance in phase-space diverge exponentially, , and thus λL > 0 is a signature of chaotic dynamics.

Formally, the Lyapunov exponent is then defined by taking the limit[36]

As the phase-space of our co-ordinate system is bounded, we evaluate Eq. (7) using the tangent-space method[35,61]. Essentially, rather than monitoring the physical separation of a pair of initially nearby trajectories, one can instead solve for the separation in tangent space, denoted by , and substitute this distance into Eq. (7). The tangent-space separation can be dynamically computed by assuming an infinitesimal initial perturbation to a reference trajectory starting at , leading to the system of equations

Here, Φ is the fundamental matrix and . The tangent-space separation with respect to the initial point in phase-space is extracted by computing with .

As we are only interested in the maximum Lyapunov exponent, it suffices to choose the initial separation along a random direction in phase-space, and we propagate Eqs. (8) and (9) for each initial condition for sufficiently large t that our estimate of λL from Eq. (7) converges.

Connection between classical and quantum Lyapunov exponents

In our discussion of the exponential growth of FOTOCs, we have argued that λQ is intimately related to the classical Lyapunov exponent λL. Specifically, we have that . Here, we further articulate this connection using a semi-classical description of the quantum dynamics, specifically by considering the evolution in the truncated Wigner approximation (TWA)[39].

First, we remind the reader that for a small perturbation δϕ, a FOTOC can be expanded to as . A simple conclusion from this expansion is that if grows exponentially we can attribute this behavior to the variance, i.e. it must be true that .

A semi-classical explanation of this exponential growth is simplified by assuming that is an operator which is linear in the classical phase-space variables . For concreteness, let us consider as in Fig. 3 of the main text, which corresponds to αR in the classical phase-space.

Next, we consider a description of the quantum dynamics within the framework of the TWA. Here, the dynamics is computed by solving the classical equations of motion (Eq. (5)) with random initial conditions sampled from the corresponding Wigner phase-space distribution of the initial state[39]. Quantum expectation values are then obtained by appropriate averaging over an ensemble of trajectories, e.g., where the overline denotes a stochastic average. The random sampling of initial conditions serves to model the quantum fluctuations of the initial state.

For a classically meaningful initial state (i.e. a product of coherent states for the phonon and spin degrees of freedom), the fluctuations in each of the phase-space variables are typically Gaussian and centered around the expectation values of the initial state. A concrete example to illustrate this is the state considered in the main text. For each trajectory, the variable (αR) (j denoting the trajectory), for example, is sampled from a Gaussian distribution with mean zero and variance 1/4. The connection between the quantum dynamics and classical chaos is made by instead considering sampling only the fluctuations δαR about a central classical trajectory, i.e. .

Solving the dynamics of the central classical trajectory and the ensemble of fluctuations is then identical to the calculation of Eqs. (8) and (9), from which the Lyapunov exponent is calculated. In particular, the connection between quantum and classical exponents is finally made clear by evaluating the quantum variance,

As δαR is evaluated directly from Eq. (9), then we expect from our previous calculations that for a generic random perturbation, sampled according to the TWA prescription, in parameter regimes where there is classical chaos. Thus, we extrapolate that the quantum variance will grow like . Inspection of this final result shows that we should expect .

Connection between FOTOCs and RE

The connection between the FOTOCs and entanglement entropy is best established by first considering the case of the spin−phonon RE . We begin by writing the purity of the reduced density matrix explicitly in terms of the elements of the density matrix,

Our insight is that, in the case of a pure global state, the summation in Eq. (12) for the purity of the reduced density matrix can be manipulated and re-expressed as

whereare the sum of the squared diagonal elements of the density matrix and the sum over the off-diagonal coherences, respectively. Thus, we seek to understand when these latter terms can be neglected and thus the purity (and associated entropy) is expressible in terms of only the .

Firstly, there is the case of a large transverse field, and an initial state which is polarized along the direction of the transverse field with vacuum occupation, i.e., . In this case, we expect the collective spin to remain strongly polarized along the field direction. If we choose the FOTOC spin rotation axis to be along that of the initial state and transverse field, , then we have that due to the absence of initial coherences between the spin sectors in this basis, and by similar reasoning . Hence, we expect Eq. (13) to simplify so that . Identical reasoning can be applied in the normal phase (B > Bc) when the phonon detuning is the largest energy scale, such that .

The second scenario is closely related to the first. Consider an initial coherent spin state polarized along an arbitrary spin direction and vacuum phonon occupation. For arbitrary transverse field strength and on sufficiently short time-scales , then the spin component of the evolved state remains largely polarized along a particular axis dictated by the initial state. Similar to the first scenario, by choosing the spin rotation of the FOTOC, , to match the polarization of the initial state, then we will have . This is justified as and vanish, as again the state at short times will not have appreciable coherences between different spin sectors in this basis.

Lastly, for a small transverse field, , and beyond short times (i.e., beyond the time-scale when the spin state is still strongly polarized and the second scenario is still valid), we expect Eq. (13) to be well approximated by for any spin rotation axis . This is because initially pure states which are sufficiently scrambled after a quench of the system parameters closely resemble so-called canonical pure thermal quantum (cTPQ) states[62] in a generic basis. For cTPQ states, the summation over off-diagonal coherences vanishes exactly for a sufficiently large system as the coherences can be considered as random variables[62]. Moreover, for a typical spin rotation axis , the cTPQ state will have a spin distribution which is largely delocalized implying that where nph is some constant which characterizes the spread of the boson number distribution. The term is then typically much smaller in magnitude when compared to the remaining terms and . This reasoning leads to . Discussion of the sensitivity of these arguments to the rotation direction can be found in Supplementary Methods.

More generally, we can extend these arguments to extract a correspondence with the Renyi entropy of a generic bipartition of the spin−phonon system. Specifically, splitting the system into a subsystem A: L spin-1/2s, and its complement A: N − L spin-1/2s and the bosonic mode. In the weak-field regime , . Here, the terms and are obtained as the Fourier amplitudes of fidelity OTOCs for generalized rotations within each subsystem. Specifically, a local rotation taken to act on the spin-12s in the A subsystem, and a joint (but uncorrelated) rotation of the spins and bosons in the complement A.

Experimental implementation

By preparing an initial spin polarized state, recent experiments[18] demonstrated it was possible to measure the many-body overlap of the final state with the initial configuration by fluorescence detection. The Dicke model, however, includes spin and phonon degrees of freedom.

While the full spin–phonon fidelity measurement has not yet been demonstrated experimentally, such measurement is possible by extending the method in[57] to a multi-qubit system. In particular, we note that this proposal is comprised of a two-step measurement, where we first measure the spin degree of freedom. The probability of all ions being in the dark state (i.e. all in the state |↓〉) can be measured with excellent fidelity and has been previously demonstrated[18]. The dark state does not scatter photons, and as such, this measurement will not change the state of the phonons. Next one can proceed to measure the phonon occupation via the protocol described in ref. [57].

Finally, as noted in the main text, we have taken into account the single-particle decoherence present in the experiment. The results presented in the main text accounted for this by approximating the effects of decoherence by an exponential decay, . We have justified this approximation by comparing to an efficient numerical solution of the full Lindblad master equation[19,63,64] for smaller system sizes (N = 10). We find that the decoherence is well-captured by the approximate model for all transverse field strengths B considered.

Thermal and diagonal ensembles

The canonical thermal ensemble, used in Fig. 4, is defined by the density matrix , which is characterized by the inverse temperature β = 1/(kBT). This inverse temperature is chosen such that energy of the ensemble is matched to that of the initial state of the dynamics, . The RE for bipartitions of the thermal ensemble is then obtained via the definition where is the reduced density matrix obtained after tracing out the phonon degree of freedom and the remaining N − L spins.

A related concept is the diagonal ensemble [1,65], which generically describes the (time-averaged) observables of a quantum system which has relaxed at long times. The ensemble is defined as the mixed state where and |E〉 are the eigenstates of the Hamiltonian with associated eigenvalue E. We use this diagonal ensemble as a comparison to the time-averaged distribution functions P(Mz) and P(n) in Fig. 4.