Identifying and harnessing dynamical phase transitions for quantum-enhanced sensing

We use the quantum Fisher information (QFI) to diagnose a dynamical phase transition (DPT) in a closed quantum system, which is usually defined in terms of non-analytic behaviour of a time-averaged order parameter. Employing the Lipkin-Meshkov-Glick model as an illustrative example, we find that the DPT correlates with a peak in the QFI that can be explained by a generic connection to an underlying excited-state quantum phase transition that also enables us to also relate the scaling of the QFI with the behaviour of the order parameter. Motivated by the QFI as a quantifier of metrologically useful correlations and entanglement, we also present a robust interferometric protocol that can enable DPTs as a platform for quantum-enhanced sensing.


I. INTRODUCTION
The isolation and control of quantum systems at the single-particle level in atomic, molecular and optical platforms has driven a surge of experimental interest in studying non-equilibrium phenomena. As a consequence, it has become clear that non-equilibrium quantum systems that feature coherence, entanglement and correlations, can be important platforms for next-generation quantum technologies [1,2].
In this manuscript, we theoretically demonstrate that the quantum Fisher information (QFI) [31], which quantifies metrologically useful entanglement and correlations in a quantum state [32][33][34], can be used to characterize DPTs via an underlying connection to excited-state quantum phase transition (EQPTs) [12,28]. Our method shares analogies with studies of the fidelity susceptibility [35][36][37] for ground-state quantum phase transitions (QPTs), but is distinguished by the addition of time as a relevant variable. We employ a numerical study of the dynamical phase diagram of the paradigmatic Lipkin-Meshkov-Glick (LMG) model to illustrate our arguments, and use a semi-analytic model to establish a quantitative connection between the scaling of the QFI and the or-der parameter. Furthermore, we demonstrate that these quantum signatures of the DPT can also be accessed through a related many-body echo combined with simple global measurements. This latter result, in particular, opens a realistic path for the harnessing of DPTs for quantum-enhanced sensing [38].

II. SIGNATURES OF DPTS IN THE QFI
To outline our arguments most generally, consider a many-body Hamiltonian describing a closed system, where [Ĥ 0 ,Ĥ 1 ] = 0 and λ is a tunable parameter. The evolution of an initial state |ψ 0 underĤ(λ) is given by |ψ(λ, t) = e −iĤ(λ)t |ψ 0 , and a time-averaged order parameterŌ = 1 T T 0 Ô (t) dt distinguishes ordered (Ō = 0) and disordered (Ō = 0) dynamical phases. A DPT is signaled by non-analytic behaviour inŌ at a critical point λ cr separating the phases [21][22][23][24]. We propose to characterize the DPT by probing how the state |ψ(λ, t) abruptly changes as the system is quenched through λ cr . This mirrors uses of the fidelity susceptibility to quantify an abrupt change in the ground-state wavefunction at an equilibrium transition [35][36][37][39][40][41]. We similarly define the QFI as the susceptibility [31,34,42], where F(λ, δλ, t) = | ψ(λ, t)|ψ(λ + δλ, t) | is the overlap between two dynamical states that differ by a perturbation δλ to the driving parameter, equivalent to a Loschmidt echo (LE) [43][44][45] F(λ, δλ, t) ≡ | ψ 0 |e iĤ(λ)t e −iĤ(λ+δλ)t |ψ 0 |. About the critical point, λ ≈ λ cr , we expect F to abruptly decrease as the dynamical states lie in distinct phases and become rapidly orthogonal, and we predict that the DPT is signaled by a corresponding sharp peak in the QFI. This signature complements the time-averaged order parameter and, similar to an equilibrium fidelity susceptibility, has the capacity to carry more information about the system as it based on a state overlap [46]. Moreover, the QFI could potentially provide a robust, agnostic probe of more sophisticated DPTs without requiring any a priori knowledge of, e.g., the appropriate order parameter for the transition. This use of the LE for DPTs is distinct from DQPTs, wherein the LE arises as a return fidelity, ψ 0 |e −iĤ(λ)t |ψ 0 , that relates to non-analyticities at a critical time rather than parameter. Moreover, our LE compares states related by a small perturbation to the Hamiltonian, whereas the survival fidelity is defined relative to a stationary state |ψ 0 and can be highly nonperturbative.

III. LMG MODEL
We demonstrate the validity of our prediction using the LMG model [47,48] as a paradigmatic example of a DPT [49]. Our choice is motivated by the collective nature of the model, which describes an ensemble of N mutually interacting spin-1/2 particles subject to transverse and longitudinal fields, as this facilitates a tractable analysis of the dynamics across a range of parameter regimes and system sizes. Moreover, the LMG model has been studied in the context of trapped ions [7,8], cavity-QED [13] and cold atoms [50][51][52][53]. It is defined by the Hamiltonian [13] whereŜ α = N i=1σ α,i /2 with α = x, y, z are collective spin operators and σ α,i are Pauli matrices for the ith spin-1/2 particle. The Hamiltonian conserves the total spin,Ŝ 2 =Ŝ 2 x +Ŝ 2 y +Ŝ 2 z , and we focus on the maximally collective sector, i.e., states with Ŝ 2 = N (N/2 + 1)/2.
The dynamical phase diagram in the classical limit (N → ∞) is shown in Fig. 1(a), for an initial state of all spins polarized along −ẑ. A pair of dynamical phases are defined in terms of a time-averaged order parameterS z and most easily described in the limit ω = 0: For Ω χ interactions force the spins to remain closely aligned to −ẑ andS z = 0, while for Ω χ the dynamics is dominated by single-particle Rabi flopping of each spin-1/2 about the +x-axis and thusS z = 0. A critical point separates the phases at Ω cr /χ = 1/2. Similar analysis holds for ω = 0, although the DPT smooths out to a crossover for ω/χ ≤ −1/8 [13,54]. The dynamical phase diagram is symmetric for Ω → −Ω in this case.
In Figs. 1(b)-(f) we use an efficient Chebyshev expansion algorithm to integrate the dynamics of a system with N = 10 3 [54, 55] and investigate the DPT using the QFI. We consider independent perturbations of the longitudinal δωŜ z [F Q,z (Ω, ω, t)] and transverse (f) Same as (e) but for fixed ω/χ = 0 and varying the tipping angle θ (relative to +z) of the initial spin state with fixed azimuthal angle φ = 0.2π. These initial states break the symmetry Ω → −Ω of the phase diagram and so we plot critical points for both positive and negative transverse field strength. Insets of (e) and (f) show scaling of Ω * with N for ω = 0 and θ = 0.9π, respectively. δΩŜ x [F Q,x (Ω, ω, t)] fields, and use an equivalent initial state |ψ 0 = |(−N/2) z whereŜ z |m z = m z |m z . Panel (b) shows typical dynamical behaviour of F Q,x /N t 2 as Ω/χ is varied and ω/χ = 0. The normalization is chosen to absorb the expected long-time growth of F Q,x ∝ t 2 . Around the critical point Ω cr /χ = 1/2 we observe a pronounced increase in the QFI in both transient and longtime (χt 1) regimes, such that F Q,x /N t 2 1. Within the transient time regime, the F Q /N t 2 shows multiple peaks for several medium values of χt that are determined numerically [54]. The QFI also distinguishes the ordered phase, F Q,x /N t 2 → 0, from the disordered phase, F Q,x /N t 2 ≈ 1. Panels (c) and (d) show our key result: The DPT as a function of Ω and ω is identified by demonstrable peaks in the QFI at long times χt 1. The change of the transition to a crossover for ω/χ < −1/8 is also indicated clearly by a broadened and diminishing peak of the QFI. Panel (e) shows that the critical value Ω * , determined from the peak positions of F Q,x or F Q,z , agrees excellently with Ω cr determined fromS z in the N → ∞ limit, up to finite-size effects (see inset).
Similarly, panel (f) demonstrates that the QFI reproduces the signature dependence of the DPT on the initial state [13]. We compare Ω * and Ω cr as a function of the initial state |ψ 0 = |θ, φ where θ is the tipping angle of the collective spin with respect toẑ on the collective Bloch sphere and we fix the azimuthal angle φ = 0.2π. For these initial states the phase diagram is no longer symmetric for Ω → −Ω, so we plot both relevant values of the critical transverse field (fixing ω/χ = 0).

IV. CRITICAL SCALING OF THE QFI
The scaling of the QFI with system size is important for identifying how any generated correlations and entanglement can be useful for quantum sensing. Specifically, the QFI provides a lower bound for the accuracy to which the driving parameter δλ can be determined, ∆λ ≥ 1/ F Q (λ, t) [31]. The standard quantum limit, e.g., the sensitivity that can be attained with quasi-classical uncorrelated states, sets a bound (∆λ 56,57]. Supralinear scaling of the QFI with N near the critical point would indicate potential uses for DPTs in quantum-enhanced sensing. We empirically fit the maximum of F Q over Ω for fixed ω/χ = 10 −4 and χt 1, i.e., F Q (Ω * ) = aN b . A finite ω is chosen to purposely break a parity symmetry of the Hamiltonian [54,58]. Figure 2(a) shows the exponent b as a function of the tipping angle θ of the initial state with fixed φ = 0. Excluding a region near the equator (θ ≈ π/2) we observe approximately constant values of b ≈ 1.5 and 1.75 for F Q,x and F Q,z , respectively, indicating the presence of metrologically useful correlations and entanglement for sub-SQL sensing, e.g., (∆λ) 2 = F Q < (∆λ) 2 SQL for large N , and which cannot be captured by mean-field theory [54].
To understand the scaling of the QFI we consider the generic HamiltonianĤ =Ĥ 0 + λĤ 1 and use an approximate analytic expression for the long-time secular contribution [34,54], Here, |n are the eigenstates ofĤ, c n = n|ψ 0 is the projection of the initial state into the eigenbasis and H nn 1 = n|Ĥ 1 |n . Equation (4) is valid when the spectrum ofĤ is non-degenerate or in the case where the degenerate eigenstates occupy different sectors of a symmetry which leavesĤ 1 invariant [54]. Numerical evaluation of Eq. (4) for the LMG model (Ĥ 1 =Ŝ z orŜ x ) agrees well with our numerical simulations of the dynamics for θ away from the equatorial plane. Note we chose a small ω/χ = 10 −4 to ensure the spectrum is non-degenerate for the case of F sec Q,z , although this is not crucial to our results [54]. Minor disagreement between F sec Q,z and F Q,z , and similarly for F sec Q,x and F Q,x near the equator, is due to corrections from transient terms ignored in Eq. (4).
Equation (4) gives a simple interpretation for the longtime QFI: A peak in the QFI is a result of enhanced fluctuations in H nn 1 attributable to either a sharp change in the properties of the eigenstates or the projection of the initial state into the eigenbasis, both of which can be correlated with the emergence of a DPT [59,60]. For many systems [12,28,61], including the LMG model [60], the DPT is triggered by the former effect due to an EQPT which leads to non-analytic features in S nn x and S nn z [58,62], as shown in Fig. 2(b) and (c) for a representative calculation with N = 1000. A sharp cusp (kink) is observed in S nn x (S nn z ) at a critical energy E cr /N = −Ω/2 for Ω/χ < 1. This is precisely the average energy of the initial state, E 0 = ψ 0 |Ĥ|ψ 0 , as Ω/χ (or also ω/χ in general) is tuned through the DPT at Ω cr /χ ≈ 1/2. Thus, the non-analytic behaviour of the DPT already corresponds closely with that of S nn z via the relation S z ≡ n |c n | 2 S nn z at long times [60]. By examining the projection of the initial state (c n ) for Ω ≈ Ω cr in Figs. 2(b) and (c), we conclude that the distribution of relevant S nn x (S nn z ) in Eq. (4) straddles the cusp (kink) and leads to the sudden increase of the QFI at the DPT. Moreover, Eq. (4) combined with the approx- [63], and ii) the energy fluctuations of an initial coherent spin state typically scale as ∆E ∼ √ N , allows us to qualitatively predict the scaling of the QFI as . From numerical diagonalization of the LMG Hamiltonian we obtain γ x 0.495 ± 0.005 and γ z 0.25 ± 0.02 which is consistent with the approximate scaling of F Q,x /t 2 ∼ N 1.5 and F Q,z /t 2 ∼ N 1. 75 . The scaling of F Q,z is intimately related to the scaling of the order parameterS z ∼ (Ω − Ω cr ) γz [care of approximation (i)] for a finite size system [64]. Thus the QFI, which is a detailed measure of how rapidly the dynamical state changes across the DPT, is intuitively governed by the sharpness of the DPT in terms of the time-averaged order parameter.
Our analysis further supports that the QFI correctly diagnoses the DPT. For N → ∞ the relative energy fluctuations ∆E/E of the initial state vanish and F sec Fig. 1(e) that show Ω * computed from the QFI approaching Ω cr from below (above)].
We note that while we understand the connection between the QFI and DPT as being facilitated by the EQPT, it should be distinguished from the latter as being a non-equilibrium result. For example, while EQPTs also feature a divergent fidelity susceptibility (and thus associated QFI), this is a property of eigenstates of an equilibrium system. In contrast, Eq. (4) would give a QFI of zero for an eigenstate of the system, as it is specifically the fluctuations of the initial state (i.e., distribution across the eigenstates) that drives the large QFI and, moreover, are inextricably related to the scaling we observe. Furthermore, directly harnessing the QFI of an EQPT for quantum sensing would be faced with the tandem challenges of controllably preparing an excited state of a many-body system and adiabatically tuning system parameters through a phase boundary. The latter inevitably becomes difficult with increasing system size as the required time scales diverge as a result of the vanishing energy gap.
We note that the above protocol assumes temporal control over both λ and δλ. This is a reasonable approach if the goal of the protocol is simply to characterize the QFI and thus the DPT. However, this assumption is not suitable from the perspective of using the DPT for metrology, where the perturbation δλ is intrinsically unknown. As a result, it is also reasonable for us to focus on the most general scenario where one does not have temporal control of δλ. Thus, we consider an alternative, but entirely equivalent, echo sequence where the perturbation is always present [ Fig. 3(a) This adapted protocol only requires that the sign of the control parameter λ andĤ 0 can be varied, while the perturbation δλ is identically present in both periods of evolution.
Our protocol is relevant for estimating the small perturbation δλ when one separately assumes that λ is well characterized and independently calibrated to suitably good precision. This is well motivated by the fact that while δλ and λ contribute to the Hamiltonian via the termĤ 1 , they may be generated by different physical effects. For example, in the case of the LMG model H 1 =Ŝ z corresponds to an energy splitting in an ensemble of two-level systems, which could be generated via distinct physical mechanisms such that there is some well controlled part λ and some unknown perturbation δλ due to stray external fields. This is no different to an assumption widely used in different metrological scenarios where one tunes a sensor to work at an optimal point [e.g., using a phase offset in an SU(2) or SU(1,1) interferometer].
In addition, one could also envision δλ as not being generated by a separate effect, but instead as some intrinsic uncertainty on top of the control parameter λ due to technical noise. In this case, the operation of "flipping the sign" of λ by turning some external control knob would be imperfect and actually realized as λ + δλ/2 → −λ + δ/2 [72], which follows our proposed echo sequence. Our sensing protocol could then be used to estimate the magnitude of this uncertainty δλ.
While F is important to compute the QFI, it is also an an optimal signal to infer δλ [66]. Nevertheless, while measuring the final state overlap can be simplified by the fact that DPTs are typically studied with simple uncorrelated initial states, technical challenges, such as the detection resolution required to adequately discriminate states, can still pose a practical hurdle for many platforms.
simple (e.g., Gaussian) initial state we expect the final state after the echo to be distinguishable by relatively simple and robust observables such as mean spin projection or occupation [70,[74][75][76]. Specifically, using the quantum Cramer-Rao bound, measurement of an observableM leads to a lower bound for the QFI, [31,65] which can be made arbitrarily tight for a judiciously chosen observable. For the LMG model, and assuming an initial state with all spins orientated along −ẑ, the final state |ψ f after the echo [see the W ψ (r, φ) [73] in Fig. 3(a)] approximates a weakly displaced coherent spin state and so is distinguishable from |ψ 0 by measurement ofM =Ŝ y for either choice of perturbation.
In Fig. 3(b) we compare (∆ω) −2 Sy to the QFI for a perturbation δω and a moderate system size N = 100 (pertinent for, e.g., trapped ion quantum simulators [8,77] and tailored to our later discussion of decoherence). We observe thatŜ y is sufficient to qualitatively replicate the transient features and long-time growth ∝ t 2 of the QFI. In panel (c) we plot the maximum of (∆λ) −2 Sy /N t 2 as a function of ω/χ and demonstrate that it qualitatively re-produces the peak in the transient maximum of F Q,z /N t 2 near ω ≈ 0, which identifies the DPT. In fact, near the DPT we find (∆Ω) −2 Sy /t 2 ∼ N 1.4 and (∆ω) −2 Sy /t 2 ∼ N 1. 78 [54], closely following the scaling of the QFI. Combined with the observation that (∆ω) 2 Sy < (∆ω) 2 SQL = 1/(N t 2 ) near the DPT, our results suggest that DPTs could be realistically harnessed for quantum-enhanced sensing by combining dynamical echoes with simple collective measurement observables.
We also probe the robustness of our results to typical sources of single-particle decoherence. Using the permutation symmetry of the LMG Hamiltonian we are able to efficiently simulate the dynamics of N = 100 qubits subject to single-particle dephasing at rate Γ [54,78]. For weak decoherence Γ/χ 10 −2 (within reach of, e.g., current state-of-the-art trapped ion quantum simulators [77]) strong signatures of the DPT remain in (∆ω) −2 Sy , even though we become limited to transient time-scales. Moreover, the sub-SQL sensitivity near the DPT remains robust in the same regimes.

VI. CONCLUSION
We have theoretically demonstrated that the QFI can be used to diagnose DPTs. While we establish a semiquantitative understanding of the QFI via an underlying connection to EQPTs, our analysis demonstrates it is an intrinsically non-equilibrium effect. Despite our focus on the LMG model our results can be widely applicable and it will be interesting to apply our analysis to a broader range of known DPTs [12,28,61,79]. Moreover, our interferometric protocol combining dynamical echoes and measurement of simple observables demonstrates that DPTs could be a promising path for sub-SQL sensing in non-equilibrium many-body systems [80][81][82][83], and one which sidesteps typical challenges, such as divergent time-scales, associated with quantum sensors based on equilibrium QPTs.

ACKNOWLEDGMENTS
We acknowledge helpful discussions with D. Barberena. Numerical calculations for this project were performed at the OU Supercomputing Center for Education & Research (OSCER) at the University of Oklahoma.

Appendix A: Quantum Fisher information
In this section we give several useful expressions for the QFI in the context of general Hamiltonian dynamics. Our discussion includes relevant details connecting the scaling of the QFI to features of the energy spectrum and we also provide an illustrative proof of the upper bound of the QFI for uncorrelated spin states.
In the limit of t → ∞, the sinc function enforces that only terms with ∆ nk = ∆ km = 0 survive in Eq. (A6), leading to Assuming the spectrum ofĤ is non-degenerate (see discussion below), Eq. (A7) can be reduced to a single sum, in which we have neglected all the sub-leading order terms due to finite time. As a result, an obvious but important condition for the validity of Eq. (A8) is thus that the coefficient of the t 2 term is non-vanishing, such that at sufficiently long times we can justifiably ignore those transient contributions of Eq. (A7). The validity of Eq. (A8) does extend to the case wherê H possesses a degenerate spectrum, if the degenerate states occupy different sectors of a symmetry that leaveŝ H 1 invariant. Considering the LMG model illustrates this statement clearly: Degenerate pairs of eigenstates occur in the case ω = 0 in the self-trapped phase due to the fact that the Hamiltonian possesses a parity symmetry generated by the spin flip operatorP = Π N i=1σ x,i . We denote the n-th pair of degenerate eigenstates by |+, n and |−, n where the "±" labels the eigenvalues ofP . For F Q,x we have that [P ,Ŝ x ] commute and so the degenerate terms in Eq. (A7) vanish and only the terms with ±, n|Ŝ x |±, n survive and lead to Eq. (A8). On the other hand, for F Q,z one should instead consider eigenstates that mix parity, (|+, n ±|−, n )/ √ 2, and for which the symmetric/anti-symmetric combination are uncoupled byŜ z and thus Eq. (A8) again applies. So that we can consider only a single basis when discussing Eq. (A8), we include a small ω/χ = 10 −4 to break the degeneracy when discussing results for ω → 0. We have confirmed using numerical calculations that the results obtained with ω/χ = 10 −4 for Eq. (A8) match excellently with those of a full numerical computation of F Q,z for arbitrarily small ω, and the special point ω = 0 does not change any of our qualitative conclusions. Panels (a) and (b) of Fig. 4 show typical time-traces of the normalized QFI F Q,x /(N t 2 ) and F Q,z /(N t 2 ) for an initial state with all spins polarized along −ẑ. In the long time limit, F Q,x /(N t 2 ) approaches a constant for all parameters, consistent with the form of Eq. (A8). On the other hand, we find that F Q,z does not demonstrate t 2 scaling for certain parameter regimes. For example, in the limit of large Ω χ, ω the Hamiltonian is dominated by the contribution of the transverse field,Ĥ ≈ −ΩŜ x , and we expect the energy eigenbasis to be close to the eigenstates ofŜ x . Hence, to leading order the diagonal matrix elements S nn z vanish for all n. As shown in Fig. 4(b), we only observe F Q,z ∝ t 2 when Ω Ω cr [red dotted and orange solid lines in Fig. 4(b)]. For Ω Ω cr [blue dashed line in Fig. 4(b)], the behaviour of F Q,z is dominated by transient corrections ignored in Eq. (A8) at the time-scales we can probe.
Consistent with this discussion, Eq. (A8) captures the dynamical phase diagram at long-times in almost perfect agreement with F Q,x (obtained via numerical calculation of the full dynamics). Conversely, while Eq. (A8) does not completely match F Q,z at long times for Ω Ω cr it nevertheless still captures the signatures of the DPT in F Q,z . In both cases, we find the scaling of the QFI with system size near the DPT, Ω ≈ Ω cr , is well captured. Specifically, by directly computing Eq. (A8) in a window of N ∈ [100, 2000] we obtain F sec Q,x ∼ N 1.5 and F sec Q,x ∼ N 1.75 which closely agree with results obtained from F Q,x and F Q,z for the same system sizes.
Insight can also be gained into the short-time behaviour of Eq. (A5). Using the Taylor expansion withĤ 1 =Ŝ x andŜ z for perturbations along Ω and ω, respectively, then the leading order behaviour of the QFI in the short-time limit is most generally Thus, for an initial state of θ = π we have F Q,x = N t 2 to leading order in time, which is consistent with the results of Fig. 5(a) at short times (marked as Region I) where F Q,x (Ω * )/N t 2 ≈ 1. In contrast, when considering F Q,x for the same initial state the relevant fluctuations vanish, [∆Ŝ z (0)] 2 = 0 and so we must retain higher order contributions in t. Specifically, the leading order behaviour of the QFI scales as t 4 and is generated by the commutator it[Ŝ z ,Ĥ]/2 = −ΩŜ y t/2 in Eq. (A9). Then, we obtain F Q,z = Ω 2 N t 4 /4, which is again in good agreement with the short-time results of Fig. 5(b) (see Region I). In particular, F Q,z (Ω * )/N t 2 collapses for the range of N considered and shows a power law dependence on t.
An intermediate regime separating the long-and shorttime limits (labelled as Region II in Fig. 5) also exists, wherein F Q,x/z (Ω * )/N t 2 shows transient oscillations that depend intimately on the structure of the spectrum and the initial state. Moreover, we identify a critical time t * x/z where the QFI F Q,x/z (Ω * )/N t 2 tends to a maximum value, corresponding to the first peak in Region II. x/z as a function of the system size N . We observe that t * x/z exhibits only a relatively weak dependence on N and that t * x,z χ ∼ O(10) in both cases for N ∈ [100, 1000].

Approximate model for scaling with system size
The scaling of the QFI can be directly traced to the emergence a non-analyticity in the energy spectrum of the LMG model. Here, we present a supporting calculation for this discussion.
Consider Eq. (A8) in the limit of large N such that one can make a continuum approximation H nn 1 → H 1 (E) for E ≡ E nn . Then, by recognizing that the QFI according to Eq. (A8) is proportional to the characteristic variance of H 1 (E), we argue that for an initial state with well defined mean-energy E and energy fluctuations ∆E E the QFI can be approximated as In the case of the LMG model, the divergence of the QFI arises due to a sharp cusp in S x (E) or kink in S z (E) at a critical energy E cr = N Ω/2 [see Figs. 2(b) and (c) in the main text]. Near the critical energy we find both observables are well described by [63] where we have used that the energy E and S x/z (E) are extensive observables and thus can be normalized to remove any dependence of A x,z , B x,z and γ x,z on system size N . Substituting Eq. (A12) into Eq. (A11), evaluating the derivative at E = E cr ± ∆E, and using that ∆E ∼ √ N for a coherent spin state, we obtain We then numerically diagonalize the Hamiltonian for a large system (N = 500) and fit S x/z (E) using Eq. (A12) near E cr to obtain γ x = 0.495±0.005 and γ z = 0.25±0.02, respectively, where the uncertainties reflect rms error in our fit. To be concrete, we fit S x/z (E) in the energy window of [E cr , E cr + n∆E] and [E cr − n∆E, E cr ] where n is a constant that we tuned through n = (1, ..., 5) to confirm the estimated scaling parameters γ x,z are stable. Example fits are shown in Fig. (S3). As commented in the main text, for the case of the LMG model the sharp features in S x,z (E) are related to a known excited-state quantum phase transition (although these are typically fundamentally distinct phenomena [28]). Thus, we highlight that in fact we expect the divergence of S z (E) near E cr is logarithmic in the thermodynamic limit, identical to the order parameter S z . However, our numerical calculations are limited to system sizes where finite size contributions will dominate and it is not possible to distinguish signatures of the logarithmic divergence.

Bounds on the QFI
The standard quantum limit for sensing can be recast in terms of the QFI as F cl Q ≤ N t 2 . This bound is con- The fits of (a) Sx(E) and (b) Sz(E) normalized by N/2 in the diagonal ensemble near the critical point. The red solid lines correspond to a fit using the formula given in Eq. (A12). The relevant fitted parameters are γx = 0.495 ± 0.005 and γz = 0.25 ± 0.02, respectively. Calculations are for N = 500, Ω/χ = 0.5, and ω/χ = 10 −4 , but we also check for robustness with N . Points within one standard deviation of the total energy for an initial state with θ = π and the same Hamiltonian parameters are used to perform the fitting (see text).
ventionally obtained by considering only quasiclassical initial states that feature no quantum correlations or entanglement and are then subject to evolution under only the driving term ofĤ LMG , e.g., ΩŜ x or ωŜ z [56,57]. Restricting to spin-1/2 systems, coherent spin states satisfy the former condition and a straightforward calculation demonstrates that for suitable choice of the tipping (θ) and azimuthal (φ) angles they saturate F cl Q ≤ N t 2 where N is the number of spin-1/2 particles.
In the main text we demonstrate that the interplay of interactions ∝Ŝ 2 z with the driving terms enables one to surpass the SQL near a DPT. Here, we emphasize that this result is a consequence of strong correlations and entanglement in the dynamically generated quantum states, rather than the nonlinearity of the dynamics at the classical level, by explicitly proving that the QFI remains bounded F Q ≤ N t 2 in the mean-field limit.
Plugging Eq. (A14) into the definition of F Q (λ, t) [ Eq. (A5)] we obtain where Var(Ô) = Ô 2 − Ô 2 and Cov(Ô,Ô ) = ÔÔ − Ô Ô correspond to the variance and covariance, respectively. In fact, Eq. (A15) is equivalent to expanding the QFI with respect to system size N since the first, the second, and the third term inside the square bracket typically scales as N , N 2 , and N 3 , respectively, which is related to the nature of effective one-, two-, and threebody interactions.
For an uncorrelated initial state, e.g., where |Ψ sp is some single-particle state, the second term i =j a α i a β j Cov(ĥ α i ,ĥ β j ) vanishes. However, the third term Cov(ĥ α i ,ĥ β jĥ γ k ) can still have non-zero contributions (e.g., for i = j, i = k, or j = k). Nevertheless, if we combine an uncorrelated initial state with the assumption that the Hamiltonian is single-body, i.e., can be decomposed asĤ ≡ j,α H α jĥ α j then only the linear terms in Eq. (A14) survive (e.g., b βγ jk (t) = 0 strictly).
In general, the coefficients and the single-particle variance are bounded by α (a α i ) 2 ≤ 1 and Var(ĥ α i ) ≤ ∆h 2 max /4 where ∆h max is the difference between the largest and smallest eigenvalues ofĥ. This leads to F Q ≤ N t 2 ∆h max [57]. For a spin-1/2 system, we haveĥ α i = σ α i /2 and thus ∆h max = 1, leading to the result F Q ≤ N t 2 .
This discussion is illustrated by using the specific example of the LMG model. We assume an initial (uncorrelated) coherent spin state of N spin-1/2 particles and consider the dynamics generated by the mean-field HamiltonianĤ with (time-dependent) coefficients (a, b, c) Rigorously,Ĥ MF is the effective Hamiltonian consistent with the equations of motion obtained fromĤ LMG and invoking a mean-field approximation.
We then compute the QFI for a generic single-body perturbation, F Q,α = 4t 2 [∆(Ŝ α ) t ] 2 . Formally, where T is the usual time-ordering operator. However, the form ofĤ MF means that Eq. (A17) can always be expressed as a simple sum over collective spin operatorŝ where A α (t), B α (t), and C α (t) are real numbers that satisfy Using Eq. (A18) we can thus express the time-average where the normalization factor N α is Here, n α is a unit vector aligned along the axis of the perturbationŜ α , and we have used the Cauchy-Scwharz inequality to obtain the second last line. Plugging Eq. (A20) into Eq. (A5) we finally obtain where the average is taken with respect to an arbitrary coherent spin state polarized in the direction of the unit vector n ψ . This result emphasizes that the nonlinear dynamics demonstrated by the classical model are not sufficient to generate the large QFI we observe in the full quantum dynamics, and instead this arises because of complex features that are generated in the quantum noise (see, e.g., Fig. 3(a) of the main text).
where J n is the nth Bessel function. The free parameters E min and E max are chosen such that the spectrum ofĤ is appropriately encompassed by the energy window [E min , E max ]. Throughout the manuscript we choose E max = N (χ + |ω|) and E min = −N (χ + |ω|).
To efficiently construct the complex Chebyshev polynomial φ n (−iĤ norm ) = (−i) n T n (Ĥ norm ) where T n (x) is the Chebyshev polynomial of the first kind, we use the recursion relation with the initial condition φ 0 = 1 and φ 1 = −iĤ norm .
The Chebyshev expansion is expected to converge exponentially with increasing N cut provided the N cut is not less than (E max − E min )t/2. To safely satisfy this requirement we also choose N cut to exceed this theoretical value by 30%. We check convergence by computing the normalization of the wavefunction | ψ(t)|ψ(t) | 2 and ensure that it deviates from unity by less than 10 −10 for all runs. Importantly, this deviation is much smaller than any perturbation to the wavefunction introduced by δΩ or δω when computing the QFI.

Open system with decoherence
The dynamics of the LMG model in the presence of single-particle decoherence can be efficiently simulated by exploiting the permutation symmetry of the model, combined with the fact that the initial states we probe are fully collective (i.e., Ŝ 2 = N 2 ( N 2 + 1) for our chosen initial states). In generality, the dynamics of the open system are described by a master equation for the density matrixρ of the spin ensemble [85],  [78,86] and citations therein.
In Fig. 7 we compare the behaviour of the timeaveraged order parameterS z and scaled inverse sensitivity (∆ω) −2 Sy /N t 2 across a range of longitudinal field strengths ω/χ and decoherence rates Γ/χ. We observe thatS z is relatively robust to Γ, as the primary effect of dephasing is to damp out oscillations in the disordered phase, consistent with recent experimental observations [13]. The sensitivity is more noticeably degraded by decoherence, particularly beyond short time scales (γt 1). Nevertheless, we observe that the transient peak of (∆ω) −2 Sy /N t 2 is preserved, albeit shifted towards earlier time scales and gradually smeared out around the transition. For large Γ/χ ∼ 0.1 we observe that the peak becomes systematically shifted away from the ideal DPT at ω/χ = 0 (consistent with the behaviour ofS z ), but it remains clearly centered near ω/χ = 0 for weaker decoherence. This last comparison is not unexpected, as the most noticable features of the DPT in the QFI arise for χt 10 and thus the naive requirement γt 1 for decoherence to be perturbative translates to the condition γ/χ 0.1 for the QFI to display robust transient signatures.
Calculations for a perturbation of the transverse field yields similar results. As shown in Fig. 8 we again observe that the final state after the echo is well distinguished by measurements ofŜ y . Panel (a) illustrates the Wigner distribution W ψ (r, φ) [73], which shows qualitatively similar features to Fig. 3 of the main text although the final state is instead typically displaced along the +ydirection. Similarly, the inverse sensitivity tracks the dynamics of the QFI [panel (b)]. However, we also note that in this case we always find the maximum (transient) sensitivity is at least as good as the SQL [panel (c)], max t [(∆Ω) −2 Sy /N t 2 ] ≥ 1, as for χt → 0 the perturbation results merely in a rotation of a simple coherent spin state, which is precisely the operational definition of the SQL. Nevertheless, a pronounced peak of max t [(∆Ω) −2 Sy /N t 2 ] still reflects the underlying DPT for weak decoherence, while the dynamical phases are still delineated by max t [(∆Ω) −2 Sy /N t 2 ] = 1 (ordered) and max t [(∆Ω) −2 Sy /N t 2 ] > 1 (disordered), respectively. A complete examination of the inverse sensitivity (∆Ω) −2 Sy for varied Ω and χt is also shown in Fig. 7.
3. Scaling of the sensitivity with system size To verify a robust correspondence between the QFI and the sensitivities obtained via the echo, (∆Ω)Ŝ y and (∆ω)Ŝ y , we compute the scaling of both quantities as a function of system size. At a fixed long time χt = 10 3 we fit the computed inverse sensitivity to the function aN b and obtain (∆Ω) −2 Sy ∼ N 1.4 and (∆ω) −2 Sy ∼ N 1.78 across a window of N ∈ [100, 2000]. These results closely follow the scaling of the QFI obtained via integrating the full quantum dynamics, F Q,x ∼ N 1.5 and F Q,z ∼ N 1.78 , extracted with the same procedure. We note that (∆ω) −2 Sy has the same system size scaling compared to F Q,z but with a smaller prefactor. As for (∆Ω) −2 Sy , both the prefactor and the system size scaling are less than those for F Q,x .

Appendix C: Classical dynamical phase diagram
The dynamical phase diagram of the LMG model can be solved analytically in the classical (N → ∞) limit (see, e.g., Refs. [13,28,49]). Briefly, the classical limit is equivalent to solving the equations of motion for expectation values Ŝ x,y,z under a mean-field approximation wherein all higher-order correlations are expressed as the product of single-body terms, e.g., Ô 1Ô2 = Ô 1 Ô 1 . Assuming an initial state where all the spins are fully polarized along an arbitrary axis, i.e., S ≡ Ŝ x , Ŝ y , Ŝ z = N 2 sin(θ)cos(φ), N 2 sin(θ)sin(φ), N 2 cos(θ) , the dynamics of the mean-field observable S z (t) = Ŝ z (t) can be reduced to an equivalent model of a classical particle in a potential, described by the differential equation: Here, the effective potential V eff is and the total energy E = − N 2 χ 2 cos 2 (θ) + Ω sin(θ) cos(φ) + ω cos(θ) , is a conserved quantity.
FIG. 9. The 2D plot of the effective potential V eff (Sz)/χ 2 N 2 as functions of (a) Ω and Sz for θ = π, φ = 0, and ω = 0, (b) ω and Sz for θ = π, φ = 0, and Ω = 0.5χ, (c) θ and Sz for φ = 0, Ω/χ = 0.5, and ω = 0, and (d) φ and Sz for θ = 0.3π, Ω/χ = 0.5, and ω = 0. Figure 9 shows various cuts of V eff (S z )/χ 2 N 2 for selected parameter combinations of θ, φ, Ω/χ, and ω/χ. For all the cases shown in Fig. 9, a transition from a double well to a single well can be seen. In the double-well regime, a local potential barrier exists between the two wells, which supports a local maximum at S * z . The relation of the total mechanical energy of the initial state in comparison to the magnitude of the potential barrier controls the dynamical phase: The ordered phase corresponds to the case where the particle is confined to a single well, whereas the disordered phase corresponds to the case where the particle has sufficient energy to traverse the barrier and oscillate between both wells. For an initial state with S z = 0 the transition between different dynamical phases can be obtained as the condition for which the classical turning point of the particle coincides with S * z , i.e., V eff (S * z ) = 0.
In general, obtaining an analytic solution for the phase boundary is non-trivial, due to the quartic nature of V eff (S z ). However, analytical expressions can be obtained in two special cases. First, for ω = 0 the local maximum occurs at S * z = 0 due to a parity symmetry of the model (i.e., the dynamics is unchanged upon the transformation S z → −S z and S x → S x ), enabling a straightforward solution of the critical transverse field [28], Ω cr = ± 0.5 cos 2 (θ) 1 ∓ sin(θ) cos(φ) .
Second, for θ = π or θ = 0 the dependence on the azimuthal angle φ is eliminated and we obtain an expression for the critical transverse field as a function of the longitudinal field [13],