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.

where 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â (â y ) is the bosonic annihilation (creation) operator of the mode, andŜ α ¼ P N j¼1σ α j 2 are collective spin operators witĥ σ α j (α = 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 S 2 ¼Ŝ 2 x þŜ 2 y þŜ 2 z , 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 B c = 4g 2 /δ. For B > B c (normal phase), the ground-state is described by spins aligned along the transverse field and a bosonic vacuum. For B < B c (superradiant phase), the ground-state is ferromagnetic, hjŜ Z ji $ N=2, and characterized by macroscopic occupation of the bosonic mode (Fig. 2a). Furthermore, in the superradiant phase (B < B c ), 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][29][30][31] at a critical energy E c = −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 E c . For E > E c the spectral statistics are similar to the Wigner−Dyson distribution P W ðsÞ ¼ πs=2expðÀπs 2 =4Þ, which in random-matrix theory describes a chaotic system. For E < E c the shape of the histograms Lyapunov exponents 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 H 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, V ¼ Ψ 0 j i Ψ 0 h j, 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 is neither Wigner−Dyson nor Poissonian P P (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.
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: x ¼ ðhŜ x i; hŜ y i; hŜ z i; α R ; α I Þ, where h i denotes the expectation values, and α R(I) is the real (imaginary) part of hâi. In the superradiant phase and for mean-field energies E > E c , two trajectories initially separated by Δxð0Þ in phase-space diverge as jΔxðtÞj $ jΔxð0Þje λ L t 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/E c (see Methods). For E < E c in the superradiant phase (B < B c ) and all energies in the normal phase (B > B c ), 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 > E c and B < B c a positive exponent is found signaling chaos. Note that the state   Fig. 3 Signatures of classical chaos in quantum FOTOCs. a Initial exponential growth of the FOTOC, ½1 À F X ðtÞ=ðδϕÞ 2 and the initial state jΨ c 0 i ¼ jðÀN=2Þ x i j0i (see Supplementary Note 1 for examples of exponential growth in other states). We assume δϕ ( 1=N such that we may equivalently use varðXÞ ' ½1 À F X ðtÞ=ðδϕÞ 2 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 t Ã $ a 0 þ logðNÞ=λ Q with a 0 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 jΨ c 0 i, the same state as (a), and bottom for jΨ y 0 i jðÀN=2Þ y i j0i, here N = 10 4 particles. In both plots we observe λ Q ' 2λ L . Error bars for λ Q are a 95% confidence interval from an exponential fitted to the numerical data. Coupling g and detuning δ are same as 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 λ Q t Ã $ logN, 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Ŵ G ¼ e iδϕĜ to be a sufficiently small perturbation δϕ ( 1 ð Þand setsV to be a projection operator onto a simple initial state |ψ 0 〉, i.e.V ¼ρð0Þ ¼ jΨ 0 ihΨ 0 j. This is because in the perturbative limit δϕ ( 1, this particular type of fidelity OTOC (FOTOC) 18,19 , F G ðt; δϕÞ hŴ y G ðtÞρð0ÞŴ G ðtÞρð0Þi (such that for a pure state F G ðtÞ jhψ 0 je iĤt e iδϕĜ e ÀiĤt jψ 0 ij 2 ) reduces to 37 1 À F G ðt; δϕÞ % δϕ 2 ðhĜ 2 ðtÞi À hĜðtÞi 2 Þ δϕ 2 var½ĜðtÞ; ð3Þ where var½ĜðtÞ 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, varðĜÞ 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][42][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 replacingV by the initial density matrix jΨ 0 ihΨ 0 j !ρ 0 , 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][45][46] .
In Fig. 3a we plot the FOTOCs of a small perturbation usinĝ G ¼X ¼ 1 2 ðâ þâ y Þ starting with jΨ 0 i ¼ jΨ c 0 i. In the superradiant phase we observe that after a short time of slow dynamics, t λ $ λ À1 Q , the FOTOCs feature an exponential growth $ e λ Q t , before saturating at t Ã $ log N (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 λ Q ' 2λ L , 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][48][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.
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, S 2 ðρ A Þ ¼ ÀlogTrðρ 2 A Þ measures the entanglement between a subsystem A and its complement A c = S − A, withρ A the reduced density matrix of A after tracing over A c . Although scrambling and entanglement buildup are closely connected, they are not the same. Nevertheless, a formal relationship between the OTOCs and S 2 ðρ A Þ 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 S 2 ðρ A Þ. 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 S 2 ðρ A Þ 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Ŝ r ðe r Á SÞ; where e r is a unit vector in the Bloch sphere, satisfyingŜ r jm 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, in such a way thatρŜ r M 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ρ The terms DŜ r ;n diag ðtÞ and CŜ r ;n off ðtÞ are explicitly detailed in the Methods, but importantly DŜ r ;n; diag ðtÞ is composed purely of the diagonal elements ofρ while CŜ r ;n off ðtÞ contains information about coherences. During unitary evolution the characteristic dephasing time of the coherences is t c $ λ À1 Q , which for scrambling systems is much faster than t Ã $ λ À1 Q log N. After t c any remaining coherences are fully randomized and destructively interfere yielding CŜ r ;n off ! 0. This feature, together with the fact that for those systems also the magnitude of DŜ r ;n diag becomes much smaller than IŜ r 0 and In 0 as the density matrix spreads out over the systems degrees of freedom, allows us to approximate  Fig. 4a we show the typical behavior of the RE, S 2 ðρ ph Þ, in the two different phases for |Ψ 0 〉. First, in the normal phase (panel (i)), B ≪ B c , 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 ≪ B c we observe a rapid growth of entanglement and saturation past the transient regime. We summarize our results in Fig. 4b 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Ŝ r at each time to minimize the coherence and diagonal terms in Eq. (4) (see Methods and Supplementary Methods).
The saturation of S 2 ðρ ph Þ for B < B c 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 ρ D with purely diagonal elements (see Methods) 1-3 . These To remove finite-size effects and residual oscillations we plot a timeaveraged value S 2 ðρ ph Þ for 4 ms ≤ t ≤ 12 ms (FOTOC quantities are averaged identically). The regular and chaotic dynamics for the initial state jΨ c 0 i are clearly delineated: S 2 ðρ ph Þ % 0 for B > B c and S 2 ðρ ph Þ > 0 for B < B c 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(M z ) 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 S 2 ðρ L A Þ (black markers) as a function of partition size L A 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 jΨ c 0 i, S therm We can also investigate the growth of entanglement on different size bipartitions for B < B c . For that we split the spin system into a subsystem of size L A ≤ N and evaluate S 2 ðρ L A Þ by computing the reduced density matrixρ L A by tracing over the bosonic degree of freedom and the remaining N − L A 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 A ≪ N (see Fig. 4d). However, for L A~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 S 2 ðρÞ ¼ 0, 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, bothV andŴ should be restricted to a partition of size L A 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 manybody echo requires the application of a spin echo π pulse along e r ¼ŷ which reverses the signs ofŜ x andŜ z 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 singleparticle 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 singleparticle decoherence is modeled by an exponential decay of the FOTOC components IĜ 0 ! IĜ 0 e ÀΓNt (see Methods). The numerical calculation indicates that even with decoherence the crossover between the two regimes at B~B c 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 cnumber expectation values, i.e.,Ŝ j ! hŜ j i for j ¼ x; y; z andâ ! hâi (where we adopt α RðIÞ as the real (imaginary) component of hâi). We thus obtain an equation of motion forx ¼ ðhŜ x i; hŜ y i; hŜ z i; α R ; α I Þ, 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 Δxð0Þ ¼ jx 1 ð0Þ Àx 2 ð0Þj diverge exponentially, jΔxðtÞj % jΔxð0Þje λ L t , 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 jΔxðtÞj of a pair of initially nearby trajectories, one can instead solve for the separation in tangent space, denoted by δxðtÞ, and substitute this distance into Eq. (7). The tangent-space separation δxðtÞ can be dynamically computed by assuming an infinitesimal initial perturbation to a reference trajectory starting atxð0Þ ¼x 0 , leading to the system of equations Here, Φ is the fundamental matrix and M ij dF i =dx j . The tangent-space separation with respect to the initial point in phase-spacexð0Þ ¼x 0 is extracted by computing δxðtÞ Φδxð0Þ with Φð0Þ ¼ I.
As we are only interested in the maximum Lyapunov exponent, it suffices to choose the initial separation δxð0Þ along a random direction in phase-space, and we propagate Eqs. (8) and (9) for each initial conditionx 0 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 λ Q ' 2λ L . 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 F G ðt; δϕÞ can be expanded to Oðδϕ 2 Þ as F G ðt; δϕÞ % 1 À δϕ 2 varðĜÞ. A simple conclusion from this expansion is that if F G ðt; δϕÞ grows exponentially we can attribute this behavior to the variance, i.e. it must be true that varðĜÞi $ e λ Q t .
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 x. For concreteness, let us considerĜ ¼X ¼ 1 2 ðâ þâ y Þ 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., hXi α R 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 jΨ c 0 i ¼ jðÀN=2Þ x i j0i considered in the main text. For each trajectory, the variable (α R ) j (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. ðα R Þ j ! α cl R þ ðδα R Þ j . 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 jδα R j $ e λ L t 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 δα 2 R À δα R 2 $ e 2λ L t . 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 S 2 ðρ ph Þ. 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 are 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 IĜ 0 . Firstly, there is the case of a large transverse field, B ) B c and an initial state which is polarized along the direction of the transverse field with vacuum occupation, i.e., jΨ 0 i ¼ jð ± N=2Þ x i j0i. 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,Ŝ r ¼Ŝ x , then we have that CŜ x ;n off ðtÞ % 0 due to the absence of initial coherences between the spin sectors in this basis, and by similar reasoning In 0 ðtÞ % DŜ x ;n diag ðtÞ. Hence, we expect Eq. (13) to simplify so that Tr½ρ 2 ph ðtÞ % IŜ x 0 ðtÞ. Identical reasoning can be applied in the normal phase (B > B c ) when the phonon detuning is the largest energy scale, such that Tr½ρ 2 ph ðtÞ % In 0 ðtÞ. 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 timescales t≲λ À1 Q , 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,Ŝ r , to match the polarization of the initial state, then we will have Tr½ρ 2 ph ðtÞ % IŜ 2 r 0 ðtÞ. This is justified as CŜ r ;n off ðtÞ and In 0 ðtÞ þ DŜ r ;n diag ðtÞ 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, B ( B c , and beyond short times t ≳ λ À1 Q (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 Tr½ρ 2 ph ðtÞ % IŜ r 0 ðtÞ þ In 0 ðtÞ for any spin rotation axisŜ r . This is because initially pure states which are sufficiently scrambled after a quench of the system parameters closely resemble socalled canonical pure thermal quantum (cTPQ) states 62 in a generic basis. For cTPQ states, the summation over off-diagonal coherences CŜ r ;n off vanishes exactly for a sufficiently large system as the coherences can be considered as random variables 62 . Moreover, for a typical spin rotation axisŜ r , the cTPQ state will have a spin distribution PðMŜ r Þ which is largely delocalized implying that DŜ r ;n diag $ 1=ðNn ph Þ where n ph is some constant which characterizes the spread of the boson number distribution. The term DŜ r ;n diag is then typically much smaller in magnitude when compared to the remaining terms IŜ r 0 and In 0 . This reasoning leads to Tr½ρ 2 ph ðtÞ % IŜ r 0 ðtÞ þ In 0 ðtÞ. 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 S into a subsystem A: L spin-1/2s, and its complement A c : N − L spin-1/2s and the bosonic mode. In the weak-field regime B ( B c , Tr½ρ 2 Here, the terms I A 0 and I A c 0 are obtained as the Fourier amplitudes of fidelity OTOCs for generalized rotations within each subsystem. Specifically, a local rotation e iϕŜ r;A taken to act on the spin-12s in the A subsystem, and a joint (but uncorrelated) rotation e iϕŜ r;Ac e iθâ yâ of the spins and bosons in the complement A c .
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 |↓〉 z ) 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, IĜ 0 ! IĜ 0 e ÀΓNt . 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ρ therm ¼ e ÀβĤ D =Tr½e ÀβĤ D , which is characterized by the inverse temperature β = 1/(k B T). This inverse temperature is chosen such that energy of the ensemble is matched to that of the initial state of the dynamics, hEi therm Tr½Ĥ Dρtherm ¼ hΨ 0 jĤ D jΨ 0 i. The RE for bipartitions of the thermal ensemble is then obtained via the definition S therm 2 Àlog Tr½ðρ therm L A Þ 2 whereρ therm L A ¼ Tr ph;NÀL A ðρ therm Þ is the reduced density matrix obtained after tracing out the phonon degree of freedom and the remaining N − L A spins.
A related concept is the diagonal ensembleρ D 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ρ D P E n jc E n j 2 jE n ihE n j; where c E n hΨ 0 jE n i and |E n 〉 are the eigenstates of the HamiltonianĤ D with associated eigenvalue E n . We use this diagonal ensemble as a comparison to the time-averaged distribution functions P(M z ) and P(n) in Fig. 4.

Data availability
The source data underlying Figs. 2-4 of the main text are provided as a source data file. Additional numerical data and computer codes used in this study are available from the corresponding author upon request.