Characterizing the dynamical phase diagram of the Dicke model via classical and quantum probes

We theoretically study the dynamical phase diagram of the Dicke model in both classical and quantum limits using large, experimentally relevant system sizes. Our analysis elucidates that the model features dynamical critical points that are distinct from previously investigated excited-state equilibrium transitions. Moreover, our numerical calculations demonstrate that mean-field features of the dynamics remain valid in the exact quantum dynamics, but we also find that in regimes where quantum effects dominate signatures of the dynamical phases and chaos can persist in purely quantum metrics such as entanglement and correlations. Our predictions can be verified in current quantum simulators of the Dicke model including arrays of trapped ions.

To date, the majority of experimental investigations in this direction have been tailored towards integrable models, featuring effective infinite-range interactions between spins, which admit analytical treatments [18][19][20][21]. Richer non-integrable models have been pursued in trapped ion systems [22], but the associated complexity of the quantum dynamics limited the theoretical analysis of the DPT to small system sizes and prevented a rigorous scaling analysis. Hence, it is highly desirable to find and study DPTs in non-integrable models featuring novel non-equilibrium phenomena that are both implementable in tunable quantum simulators and theoretically tractable under controllable approximations.
We advance this direction by studying a DPT in the iconic Dicke model [23][24][25], which describes the collective coupling of many spins to a single harmonic oscillator. The model is attractive as it features an array of phenomena, such as non-integrability [26][27][28], chaos [29][30][31] and equilibrium quantum phase transitions (QPTs) in both ground and excited states [32][33][34][35][36][37], but involves only a pair of disparate degrees of freedom such that the model remains amenable to analytic and numerical treatments. Moreover, the model is already studied in state-of-the-art AMO quantum simulators, including trapped-ion arrays [38,39] and cavity-QED [40][41][42][43]. Here, we investigate the DPT in analytically tractable spin-and boson-dominated limits, as well as a non-integrable regime where near-resonant coupling of spin and bosons leads to chaotic dynamical phases also seen in other non-integrable systems [7,8]. Following recent works linking DPTs to coexisting excited-state QPTs (EQPTs) [20,36,44], we find different dynamical critical points in the spin and boson-dominated regimes that reflect the presence of distinct EQPTs in these limits. By studying large, experimentally relevant, system sizes using efficient numerical methods we are able to demonstrate that mean-field features of the dynamics remain valid in the exact quantum dynamics. Conversely, in regimes where quantum effects dominate we find that signatures of the DPT and chaos persist in purely quantum metrics such as entanglement and correlations. Model: The Dicke Hamiltonian for N spin-1/2 particles is given by [30], Here,â (â † ) is the bosonic annihilation (creation) operator of an oscillator with frequency δ,Ŝ α = 1/2 N j=1σ α j are collective spin operators for α = x, y, z andσ α j Pauli matrices of the jth spin. The spins are subject to a transverse field of strength Ω and the spin-boson coupling is characterized by g. The Hamiltonian exhibits a Z 2 parity symmetry associated with the operator Π = e iπ(Ŝx+â †â +N/2) such that [Ĥ D ,Π] = 0. The equilibrium phase diagram of the Dicke model features a ground-state QPT at a critical coupling g QPT = √ δΩ/2 [32,33], which delineates a superadiant (g g QPT ) phase, with degenerate ground-states associated with different parity sectors that in the limit Ω → 0 take the form |ψ s = 1 √ 2 [|−(N/2) z ⊗|α s ±|(N/2) z ⊗|−α s ], and a normal (g g QPT ) phase, with |ψ n ≈ |(N/2) x ⊗ |0 in the limit of large Ω [38]. Here, we have defined collective spin states viaŜ x,y,z |m x,y,z = m x,y,z |m x,y,z and |±α s is the bosonic coherent state for α s = g √ N /δ. Additionally, an excited-state QPT (EQPT) exists in the spectrum of the superradiant phase, g > g QPT , defined by a critical energy E c = −ΩN/2 [34,36,37]. The EQPT is signaled by a non-analyticity in the density of states [35] and the existence of pairs of degenerate eigenstates with different parity for E < E c . Dynamical phase diagram: We study the dynamical phase diagram that arises after a quench of the transverse field. Concretely, the system is initialized in the ground-state ofĤ D at fixed g and δ with Ω = 0, such that |ψ(0) = |(−N/2) z ⊗ |α s , and the transverse field is then quenched to a final value Ω = 0. To garner insight into the dynamics we first study the classical model (Figs. 1-3) before probing the role of quantum fluctuations. The classical limit of the Dicke model is equivalent to solving the Heisenberg equations of motion for operators under a mean-field approximation, wherein expectation values are factorized according to . For brevity we adopt the notation O ≡ Ô (t) herein.
It is convenient to study dynamics in terms of two variables:g = 2g/ √ δΩ ≡ g/g QPT and δ/Ω. The former characterizes the effective strength of the spin-boson interaction relative to the single-particle terms, while the latter describes the relative energy scales of spin and bosonic excitations and loosely expresses the relative importance of each degree of freedom to the dynamics. In the limit δ/Ω 1 andg ∼ 1 the dynamics are equivalent to a spin model described by the Lipkin-Meshkov-Glick (LMG) Hamiltonian,Ĥ eff = (χ/N )Ŝ 2 z + ΩŜ x with the boson-mediated interaction characterized by χ ≡ 4g 2 /δ [17,46]. For δ/Ω 1 andg ∼ 1 the spins are instead slaved to the dynamics of the dominant bosons [45]. However, in this limit a simple boson-only Hamiltonian is not generally applicable [47]. Lastly, for δ/Ω ∼ 1 the dynamics is complicated and involves both degrees of freedom. Herein, we refer to the integrable limits δ/Ω 1 and δ/Ω 1 as the spin-dominated (SDR) and boson-dominated (BDR) regimes respectively, and the non-integrable case of δ/Ω ∼ 1 as the resonant regime (RR).
The typical dependence of the dynamics ong is shown in the time-traces of Fig. 1 for SDR and BDR. In both limits we identify that the dynamics can be characterized as either trapped or untrapped. The former occurs when the spin-boson interaction dominates the Hamiltonian and leads to a locking of the spins and bosons close to their initial configuration. In the SDR this has been interpreted as a self-generated detuning ∝ Ŝ z Ŝ z that locks out rotations due to the transverse field [21]. Conversely, the untrapped dynamics are characterized by large coherent oscillations in S z and X = 1 2 â +â † dominated by the single-particle terms of the Hamiltonian.
We classify the dynamical phases and identify a DPT using a pair of interchangeable time-averaged order  √ 2 in the spin-dominated regime and ii)g b DPT ≈ 3 1/4 in the boson-dominated regime. The dynamical phases in the LMG model have recently been observed in a cavity-QED quantum simulator [21].
In the non-integrable RR, δ/Ω ∼ 1, the dynamical phase diagram is more complex. The main panel of Fig. 1 illustrates that the time-averaged order parameter becomes noisy in the untrapped phase. In this regime typical time-traces [ Fig. 2(a)] feature erratic oscillations in both spin and boson observables and there exist short periods where the system abruptly becomes re-trapped. This behaviour signals a chaotic dynamical phase [7,8] that arises due to the known non-integrability of the Dicke model forg > 1 and δ ∼ Ω [27,29,30,32].
We support this by computing the Lyapunov exponent λ L [30] in Fig. 2 (see also Ref. [45]), as a function ofg and δ/Ω for the initial condition corresponding to |ψ(0) . Chaos, and thus non-integrability, is signalled by λ L > 0 [48] in regions of parameter space that qualitatively coincide with the noisy behaviour of S z (Fig. 1). We also find evidence of small regions where the chaotic dynamical phases exist but λ L → 0 (within numerical error). This is not contradictory as the underpinning requirement for the so-called chaotic dynamical phase to exist is in fact the non-integrability of the Dicke model due to the coupling of the spin and boson degrees of freedom.
The dynamical phase diagram of SDR and BDR is captured by an effective model involving only the dominant degree of freedom [45, [49][50][51][52][53]. Specifically, the mean-field dynamics are reduced to an equivalent picture of a classical particle with co-ordinate ξ = S z (spin-dominated) or ξ = X (boson-dominated) confined within a 1D potential V (ξ) that depends only ong, δ/Ω and the initial state. Near the DPT V (ξ) is well approximated by a double well potential [ Fig. 3(a)].
The initial position and velocity of the particle, ξ(0) andξ(0), and the relative height of the central maximum, V (0), characterizes the dynamics of the system [8]. For g >g DPT and the particle initially located in one well, ξ(0) = 0, the particle has insufficient mechanical energy to overcome the barrier and remains confined. On the other hand, forg <g DPT the particle has sufficient energy to pass over the barrier and traverses freely between both wells. The former scenario describes trapped dynamics, ξ = 0, whilst the latter describes the untrapped phase, ξ = 0. The critical point is defined as the condition for which the particle first surpasses the central barrier and we findg s DPT = √ 2 andg b DPT = 3 1/4 , respectively, in agreement with Fig. 1.
The DPT in the SDR coincides with the well known EQPT of the Dicke Hamiltonian [34,35,37], i.e., the energy of the initial state matches the EQPT critical en- In the integrable limits the DPT can be described as a particle (green markers) in a 1D potential with co-ordinate ξ = Sz (δ/Ω 1) or ξ = X (δ/Ω 1). At smallg gDP T (untrapped phase, left) the particle has sufficient energy to traverse both wells of the potential, whereas for largeg gDP T (trapped phase, right) the particle remains energetically confined to a single well. (b) Typical trajectories in phase-space for trapped (red) and untrapped (blue) phases at δ/Ω = (0.1, 1, 10) (top to bottom). For the spin phase-space we use co-ordinates r = 1 + 2Sx/N and ϕ = arctan(Sy/Sz) (c) Dynamical phase diagram as a function of initial state. At small detuning δ/Ω = 0.1 we vary the boson amplitude α (top panel), and at large detuning δ/Ω = 4 we vary the tipping angle θ0 relative to the south pole of the Bloch sphere (Sz = −N/2) (bottom panel).
ergy, E 0 ≡ ψ 0 |Ĥ|ψ 0 |gs DPT = E c . The EQPT features a non-analyticity in eigenstate observables at E c that can be intimately related to the behaviour ofS z at the DPT and is related to a saddle point that emerges in the classical phase-space [35,37]. In the BDR the Dicke model has been predicted to feature a set of EQPTs at energies {E c } that arise in distinct sectors of the energy spectrum labelled by values of an emergent conserved quantity, which restores the integrability of the model [52,53]. In the large N limit our state occupies a single sector of the spectrum and atg b DPT <g s DPT the energy of the initial state E 0 ≡ ψ 0 |Ĥ|ψ 0 |gb DPT matches a different EQPT critical energy E , which is consistent with the observed shift in the DPT critical point.
The potential model is further evidenced by inspecting typical trajectories of trapped and untrapped cases in the classical phase-space, as shown in Fig. 3(b) for the three regimes. In BDR and SDR we observe well defined orbits in the dominant degree of freedom that are centered around fixed point(s) in phase-space corresponding to the minima of the potential well(s) [45]. The trajectories of the complementary slaved observables show similar behaviour, but tend to densely fill out the available phase-space due to fast micromotion on top of the slower orbits arising from the enslavement to the dominant degree of freedom (see also Fig. 1). In the resonant regime the potential description breaks down as no single degree of freedom dominates. Consequently, we observe a lack of clear orbits in phase-space, although the erratic trajectories do fill out distinct volumes of phase-space in each dynamical phase.
DPTs in isolated systems are intrinsically dependent on initial conditions [21]. We demonstrate this by probing initial states of the form |ψ(0) = |(−N/2) n ⊗ |α where α ∈ R and we defineŜ n |m n = m n |m n for S n =Ŝ · n with n = (0, sin(θ 0 ), cos(θ 0 )) defined by the tipping angle θ 0 from the south pole of the collective Bloch sphere. In Fig. 3(c) we probe the dependence of g DPT on the initial amplitude α in the BDR, while in the SDR we vary the initial tipping angle θ 0 . We observe good agreement with analytic predictions forg DPT based on the potential model [45]. Small deviations are observed for θ ≈ ±π/2 and α/α s → 0, which correspond to initializing a particle near the central maximum of the potential such that even a small off-resonant exchange of energy with the complementary degree of freedom is enough for the particle to overcome the barrier. Quantum dynamics: Using an efficient exact diagonalization method [45] we are able to simulate the quantum dynamics and observe mean-field signatures of the DPT that survive quantum noise. For example, in Fig. 4(a) we probe a typical time-trace of Ŝ z in the resonant regime [same as Fig. 2(a)] for experimentally realistic N = 40, 200, 600 and find that it is possible to gradually observe the signature re-trapping dynamics with increasing N , before quantum fluctuations dephase Ŝ z → 0.
Conversely, pure quantum measures such as entanglement and quantum correlations also show clear signatures of the DPT, chaos and non-integrability. Figure  4(b) shows that the critical point of the DPT is signaled in the pronounced build-up of spin-boson entanglement in the untrapped phase, quantified by the entanglement entropy S vN = −Tr[ρ s log(ρ s )] whereρ s is the reduced density matrix of the spins. For δ/Ω 1 the timeaveraged entropy dtS vN (t) demarcates the trapped and untrapped phases consistent with the order parameter S z , with S vNtrapped S vNuntrapped . The entanglement decreases with increasing δ/Ω as the bosons become eliminated, although we expect that, e.g., bipartite entanglement between the spins will retain indications of the DPT. Panel (c) also demonstrates that the temporal fluctuations of the entanglement vary with the underlying integrability of the Hamiltonian [54,55]. Comparing to the Lyapunov exponent in Fig. 2, we see a correlation between regions of chaotic (non-integrable) dynamics and suppressed temporal fluctuations of S vN (t) is obtained for N = 600. Time averages in (b) and (c) are computed over a window t ∈ [t0, T + t0] for gt0, gT = 150, so as to reduce transient effects in steady state estimates. Gray regions correspond to excluded data.
(relative to S vN ).
Panel (d) evidences that the critical region of the DPT can be diagnosed by the rapid growth of quantum fluctuations. By empirically fitting (∆Ŝ z ) 2 = Ŝ 2 z − Ŝ z 2 ∝ (1−e −γt ) we find that the rate at which quantum fluctuations buildup, characterized by γ, is largest forg ∼g DPT . This is consistent with the effective potential description, as near the critical coupling the mean-field trajectory spends an increasing amount of time probing the central barrier (e.g., unstable fixed point) [8,45,56,57], motivating the expectation that quantum fluctuations will grow exponentially. Nevertheless, we also note that a similarly rapid growth of fluctuations is observed away from the DPT for 1 δ/Ω 2 as a consequence of classical chaos [30], indicating that it is important to carefully differentiate the effects of chaos from the DPT in quantum dynamics. Experimental realization: The full dynamical phase diagram of the Dicke model can be studied in a range of current state-of-the-art AMO quantum simulators, but most readily in arrays of trapped ions. In particular, the Dicke model was recently realized in a 2D Penning trap configuration [38,39] where a spin-1/2 is encoded in two internal states of each ion and the collective center-ofmass motional mode of the ion crystal realizes the bosonic degree of freedom. A pair of lasers generates an optical dipole force that couples the internal state of each ion to the motional degree of freedom [58]. A transverse field is generated by applying microwaves with tunable strength Ω that coherently drive the two internal levels and the spin-boson coupling g can be in principle controlled via the applied laser power. Moreover, δ can be varied by controlling the detuning of the lasers relative to the frequency of the targeted center-of-mass mode [58]. Considering the simulator reported in Refs. [38,39], we predict it should be possible to probe regimes 0.1 δ/Ω 10 forg ∼ 1 on timescales that are fast compared to relevant sources of single-particle decoherence of the spins [45,46,59]. In contrast to related studies of the Dicke model in optical cavities, the trapped ion simulator has little damping of the bosonic mode and thus to an excellent approximation can be neglected. Further, ion crystals of N ≈ 200 [59] are possible in the 2D geometry, which is sufficient to access the signatures of the DPT as in Figs. 1 and 4 [45]. Looking ahead, implementing the Dicke model in 3D ion crystals [60] would open a path to even larger system sizes that are beyond the capability of numerical methods.

Conclusion:
We have studied a DPT in the Dicke model and unique features arising from non-integrable and chaotic regimes of parameter space. Our numerical study indicates signatures of the DPT survive in the quantum dynamics and are accessible in current state-of-the-art AMO quantum simulators based on, e.g., trapped-ion arrays. This can motivate future investigations connecting to quantum phenomena such as, e.g., information scrambling [30,61]. Additionally, studying entanglement dynamics associated with DPTs in collective systems will open new directions for the generation of metrologically useful states for quantum-enhanced sensors [62] or frequency and time standards [63].

ANALYTIC TREATMENT OF DPT IN INTEGRABLE REGIMES
In this section we present a detailed outline of the classical particle in a potential description of the DPT in the integrable spin-and boson-dominated limits. We derive each case separately, drawing from the above mean-field EOM [Eq. (S2)], and present the dependence of the critical point of the DPT on the initial state.
Spin-dominated regime: δ/Ω 1 First, we address the limit where the spins are the dominant degree of freedom. To derive an effective model of the dynamics we take the limit δ/Ω → ∞ but keeping Ω fixed and formally eliminate the bosonic degree of freedom so we are left with a simpler pure spin description. This is achieved by observing from Eq. (S2) that the bosons will evolve on a much faster timescale (1/δ) than the spins (1/Ω). Thus, the dynamics of the bosons is composed of two different contributions: a slow drift created by the slow time evolution of s z plus very fast oscillations about it. The slow dynamics is characterized by the locally averaged valueβ, which is obtained by settingβ = 0 in Eq. (S2): As a first approximation, the spin degree of freedom will only be sensitive to the dynamics ofβ av , such that the spins evolve under a closed set of equations of motion (EOM): FIG. S1. Effective potential for the spin degree of freedom when δ Ω and θ0 = 0 (left) or θ0 = π/4 (right). The "particle" begins at sz = − cos(θ0) and then moves to the right. Wheng <g s DPT (0) (blue) the central barrier is below the initial "energy" and so sz oscillates between − cos(θ0) and cos(θ0). Wheng >g s DPT (blue) the particle cannot overcome the central barrier and so sz < 0 for all times. As illustrated by the two contrasting panels, the exact form of the potential depends on the initial state of the system (parameterized by θ0 in this instance), and in the spin-dominated regime this solely introduces the dependence of the critical couplingg s DPT (θ0) on the initial state.
These equations can equivalently be derived from a classical Hamiltonian with the standard angular momentum Poisson brackets: The energy h is a conserved quantity defined with respect to the initial condition of the spins, [s x (0), s y (0), s z (0)] = [sin θ 0 , 0, − cos θ 0 ]: Similarly, the total spin length, s 2 x + s 2 y + s 2 z = 1, is also conserved under the dynamics. Together, this allows us to eliminate s x (using conservation of h) and s y (using the EOM), to obtain a single equation for s z : This differential equation describes the motion of a classical particle in a one-dimensional potential V (s z ). The dynamics of the particle is dictated by the form of the potential and the initial condition: The particle is initialized at s z (0) = − cos θ 0 with potential energy V (s z (0)) = 0 and kinetic energyṡ z (0) 2 /2 = 0. Forg → 0, the particle can freely move from s z = − cos θ 0 to s z = cos θ 0 and back. This describes the untrapped regime. Asg is increased, V (s z ) develops a maximum at s z = 0. The height of this central barrier increases until it matches the initial mechanical energy of the particle for a critical couplingg Forg ≥g s DPT (θ 0 ) the particle is confined to s z ≤ 0 for all times (see Fig. S1). This is the trapped regime. In the case θ 0 = 0, which is the primary focus of the main text,g s DPT = √ 2.
Boson-dominated regime: δ/Ω 1 Here, we address the limit where the bosons are the dominant degree of freedom. We direct the interested reader to Refs. [49,50,52,53] for related approaches in the literature. To derive an effective model of the dynamics we take the limit δ/Ω → 0 keeping δ fixed. Complementary to the prior analysis, the spin degree of freedom now has a very short evolution timescale (1/Ω) as compared to the bosons (1/δ). Hence, from the perspective of the spins, theβ term in their equations of motion [see Eq. (S2)] is static. Hence, the spin dynamics can be understood as a fast precession with respect to a "static"β-dependent axis. Furthermore, the slow time evolution ofβ induces adiabatic evolution of the spins on top of their fast precession, with instantaneous axis of rotation given by whereβ R = (β +β * )/2 is the real part ofβ and n t owes its time dependence toβ R . The drastic difference in time-scales between the spin and bosonic degrees of freedom means that the bosons are only sensitive to the local average of the spin vector (e.g., we ignore the rapid oscillations on 1/Ω timescales) s: with s 0 ≡ s(0) the initial condition. For simplicity of analysis in the following, we will take θ 0 = 0 (relevant for the main text). Then, the boson equation of motion becomeṡ which can be recast in terms of an effective potential by eliminating the imaginary part ofβ: The analysis of this differential equation follows identically to that of the spin-dominated regime. Forg → 0 the particle moves freely between ±β 0 , which corresponds to the untrapped regime. Asg increases, a maximum develops atβ R = 0 and eventually acts as a barrier that keeps β 0 > 0 for all times. This occurs at the critical coupling (S13) For β 0 = 1, we recoverg DPT = 3 1/4 as per the result quoted in the main text. An important difference with the spin potential is that in this case the form of the effective boson potential depends on the initial state of the spins as well as the initial state of the bosons, whereas the spin potential is insensitive to the initial state of the bosons. While our result here, Eq. (S13), specializes to θ 0 = 0, we emphasize that the calculation is easily generalized for arbitrary θ 0 .

COMPARISON TO EQPT
The Dicke model has been intensely studied with regards to the presence of a transition in the excited-state spectrum of the Hamiltonian (EQPT). Not only was the nature of the EQPT subject to initial debate [34][35][36], although it is now accepted as a true quantum phase transition, but also there was initially speculation as to whether the EQPT might be connected to various non-equilibrium features of the model such as the emergence of chaos in the energy spectrum [29,34].
In general, the Dicke model features an EQPT forg > 1 (e.g., when the coupling is large enough that the groundstate corresponds to the superradiant phase) at a critical energy E c = −ΩN/2. It is signaled by a non-analyticity in the density of states, in particular a logarithmic divergence in the first derivative [35], as well as singularities in Ŝ x and â †â computed with respect to the energy eigenstates [34]. The EQPT can be explicitly connected to the classical model in the large-N limit by noting that it arises due to the emergence of a saddle point in the classical phase-space [35,37].
The DPT we study is quantitatively connected to this EQPT in the spin-dominated regime: We find that the energy of the initial state, E 0 = ψ(0)|Ĥ D |ψ(0) = −(g 2 /2)ΩN/2, coincides with the critical energy E c = −ΩN/2 at g s DPT = √ 2. Indeed, going a step further, it is straightforward to prove that the DPT in the spin-dominated regime always precisely occurs when E 0 = E c for any initial state (e.g, arbitrary θ 0 , as covered in the previous treatment of the DPT by the particle in a potential model). This close correspondence is consistent with prior work in the literature [36], which suggested that EQPTs might be accessible through a type of DPT driven by the breaking of the parity symmetry in the energy spectrum.
In the boson-dominated regime the situation is made complex by the fact the model becomes approximately integrable. Specifically, previous work [52,53] has demonstrated that the Dicke model in the limit δ/Ω → 0 features an emergent conserved quantity,Ŝ x =Ŝ which splits the energy spectrum into distinct sectors. Within each sector, labeled by the value m x of the conserved quantity, the model features an EQPT at the critical energy E b c (m x ) = −Ωm x . In the limit of large N , our initial state corresponds to the sector with and we find that the energy of our initial state at the BDR DPT then precisely coincides with the EQPT,

RELATION BETWEEN LYAPUNOV EXPONENT AND CHAOTIC DYNAMICAL PHASE
In the main text we highlight that the chaotic dynamical phase can be correlated closely with the known chaos of the Dicke model, as quantified by the Lyapunov exponent (see Figs. 1 and 4, respectively). Figures S2(a) and (b) quantify this relation more closely. Particularly, in panel (b) we overlay the dynamical phase diagram with a contour line corresponding to λ L = 0.01 in (a), approximately indicating that the interior region of the greend boundary corresponds to chaotic dynamics. We observe that this reasonably aligns with the noisy (e.g., speckled) patches of the dynamical phase diagram.
We briefly note that disagreement between the region of non-zero Lyapunov exponent and the order parameter noise, such as, e.g., for (δ/Ω,g) ≈ (1.25, 2.5) arises due to the fact that chaos is not a one-to-one indicator of nonintegrability. In turn, so-called 'chaotic' dynamical phases only require the underlying model to be non-integrable [7,8], and so the lack of chaos for a given set of parameters does not exclude the existence of 'chaotic' dynamical phases.

NUMERICAL SIMULATION OF QUANTUM DYNAMICS
We efficiently simulate the quantum dynamics of the Dicke model by solving the time-dependent Schrödinger equation using a Krylov-subspace projection method [64]. Using the symmetries of the system, we are able to make our numerical calculation tractable by: i) representing the spin degree of freedom in the fully symmetric subspace of Dicke states |S, m z with S = N/2, defined such thatŜ 2 |S, m z = S(S + 1)|S, m z forŜ 2 = α=x,y,zŜ 2 α andŜ z |S, m z = m z |S, m z , and ii) representing the boson degree of freedom in a truncated Fock basis |n where n ∈ [0, n max ] for some maximum occupation n max . While it is straightforward to simulate the evolution using a fixed n max , this may require extensive trial runs and benchmarking to ensure convergence of all relevant observables.
To avoid this issue, and to efficiently simulate large system sizes to long times over a wide range of parameter space, we utilize a bosonic Hilbert space that dynamically grows or shrinks throughout the evolution. Initially, we set n max by bounding the error in the initial wavefunction normalization 1 − | ψ(0)|ψ(0) | 2 ≤ = 10 −6 ; by checking for convergence as is varied, we have verified that this bound is sufficient to generate relative errors in the spin-boson entanglement entropy dynamics on the order of 10 −5 at long times, and even smaller errors for relevant observables. Throughout the evolution, we periodically check (on a time interval τ ) the normalization of our wavefunction when projected onto states with boson occupations of n max and n max − 1. If the error introduced by the projection onto either of these sets of states exceeds at any point in this interval, we increase n max by a fixed amount ∆n and recompute the evolution over this interval, repeating until the error tolerance is satisfied. In general, we find that the short-time dynamics requires a much larger value of n max than the initial state; in contrast, however, the late-time dynamics typically require a much smaller value of n max . To speed up the computation of the late-time dynamics, we also allow for a decrease in the value of n max , likewise checking the projection onto all boson occupations between n max − ∆n and n max . If the projection onto this set of states is less than at any point over some interval τ , we reduce n max by ∆n, and truncate away the corresponding states from our Hilbert space.
For all simulations in Fig. 4(b)-(d), we use N = 600 and evolve to a fixed time gt = 300, initializing the system in the state |(−N/2) z ⊗ |α s with α s = (g/δ) √ N as described in the main text. In Fig. S3, we show typical dynamics of the entanglement entropy, S vN = −Tr[ρ s log(ρ s )] whereρ s is the reduced density matrix of the spins, magnetization Ŝ z , and variance (∆Ŝ z ) 2 = Ŝ 2 z − Ŝ z 2 for various parameter regimes. To produce Fig. 4 of the main manuscript, we simulate the Dicke model for a range of 51 values ofg and 39 values of δ/Ω, equally-spaced between the bounds shown on the plot. For a small set of parameters, typically with δ/Ω 0.3 and largeg 1.4 (corresponding to large initial boson occupation), the required time and memory resources for simulation are much greater than other regions of parameter space. Thus, they are thus excluded from our results and indicated by the gray regions in Fig. S3.
For the entanglement entropy and spin variance in Fig. S3, we also display fits to an empirical function A(1 − e −γt ) with free fitting parameters A (steady state value) and γ, shown as dashed lines in Fig. S3. We generally find that the entanglement entropy S vN (t) is consistent with such a functional form, with the exception of transient oscillations on fast time-scales. For some parameters we find that the early time dynamics of the spin variance (∆Ŝ z ) 2 is captured more accurately by a translated logistic function, a/(1 + e −c(t−t0) ) − b, than an exponential function; here, a provides the overall scale of the function to be consistent with the observed steady state, c is the logistic growth rate, t 0 is the function midpoint, and b is an offset to allow for the fact that (∆Ŝ z ) 2 = 0 at t = 0. Fitting to this function provides an alternative estimate for the relaxation timescale of the dynamics via γ logistic = c + 1/t 0 , where we account for both the logistic growth timescale as well as the delay in the onset of this growth.
In Fig. S4, we plot various notions of the relaxation timescales for the dynamics, including the value of γ obtained by an exponential fit to the spin variance, which is plotted in Fig. 4(d) of the main text. Despite the complex nature of the dynamics over the wide range of behaviors we observe, we find that these timescales -obtained by fits to relatively simple functions -are generally consistent with each other, and thus serve as reliably proxies in systematically comparing the dynamics over the considered parameter space.
Lastly, in Fig. S5, we also provide results for N = 200 and N = 40, analogous to those for N = 600 shown in Fig. 4 of the main manuscript, with the addition of the steady-state magnetization Ŝ z . For N = 40, the timescale γ is extracted via an empirical fit of the entanglement entropy S vN (t) ∝ (1 − e −γt ), as opposed to the spin fluctuations (∆Ŝ z ) 2 , which is used for N = 200 and N = 600; the spin fluctuations for smaller systems exhibit noisy oscillations that prevent a consistent estimate of a relaxation timescale. While signatures of the DPT are still evident in the single-body observables for systems as small as N = 40, features of both the DPT and classical chaos are overshadowed by the associated quantum noise in this system. However, for N = 200, which is within the capabilities of current trapped-ion platforms, signatures of classical chaos and the DPT feature much more prominently.

EXPERIMENTAL REALIZATION IN A TRAPPED ION QUANTUM SIMULATOR
In the main text we discuss how the dynamical phases we report can be realized in an implemenation of the Dicke model in a 2D trapped ion array formed in a Penning trap experiment [38,39]. Here, we support these comments with a brief comparison of the achievable coherent time-scales relative to decoherence and technical limitations.
First, we establish the time-scales of the dynamics we are most interested in. To approach this problem concretely, we will focus on the limiting cases of the boson-and spin-dominated for which we can identify the relevant energyscales from our prior analytic treatment in the mean-field limit. We expect quantum corrections to predominantly lead to, e.g., damping of oscillations due to the quantum noise. The resonant regime interpolates between the integrable spin and boson regimes, and thus should be suitably estimated from our results here.
The key requirement of the experiment will be to distinguish the untrapped from trapped phase. In current   [38]. Solid lines indicate δ/Γ el for boson-dominated regime (δ/Ω 1), while dotted lines indicate Ω/Γ el for spin-dominated regime (δ/Ω 1). Red shaded areas are relevant for comparison with δ/Γ el : Upper red shaded area indicates detunings that may introduce couplings to other normal modes of the crystal, δ/(2π) > 10 kHz, while lower shaded area indicates detunings below the limit imposed by COM fluctuations δ/(2π) < 50 Hz.