Signatures of quantum phase transitions after quenches in quantum chaotic one-dimensional systems

Quantum phase transitions are central to our understanding of why matter at very low temperatures can exhibit starkly different properties upon small changes of microscopic parameters. Accurately locating those transitions is challenging experimentally and theoretically. Here we show that the antithetic strategy of forcing systems out of equilibrium via sudden quenches provides a route to locate quantum phase transitions. Specifically, we show that such transitions imprint distinctive features in the intermediate-time dynamics, and results after equilibration, of local observables in quantum-chaotic spin chains. Furthermore, we show that the effective temperature in the expected thermal-like states after equilibration can exhibit minima in the vicinity of the quantum critical points. We discuss how to test our results in experiments with Rydberg atoms, and explore nonequilibrium signatures of quantum critical points in models with topological transitions.


I. INTRODUCTION
Quantum phase transitions are key to our perception of quantum matter across fields in physics, from quarkgluon plasma and neutron stars to quantum magnets and high-temperature superconductors 1,2 . At those transitions different quantities in completely different systems can exhibit universal behavior. This is something that we understand thanks to the development of the renormalization group theory. Among the challenges that remain for each specific system is to (if possible) find experimentally where quantum phase transitions occur, as well as theoretically predict their locations using simplified models. Quantum simulators promise to overcome the latter challenge by experimental means, as they provide pristine and controllable realizations of theoretical models 3,4 .
Quantum simulators also provide access to real-time dynamics. This is something that can be used to explore unique aspects of crossing a quantum phase transition in real time. For example, recently a Rydberg-atom quantum simulator was used to probe the Kibble-Zurek mechanism of universal defect production for slow parameter sweeps 5 . On the theoretical side, recent works have provided evidence that nonequilibrium quantum evolution can be used to probe quantum phase transitions in integrable systems 6,7 , in prethermal states for models close to integrability 8 , or through out-of-time-order correlators 9 . However, identifying real-time signatures of quantum phase transitions in generic (quantum-chaotic) many-body systems has remained a challenge.
In this work we show that generic quantum matter can exhibit dynamical signatures of quantum phase transitions by the antithetic strategy of forcing these systems out of equilibrium and therefore beyond the ground-state manifold. We find that the intermediate-time dynam-ics of local observables and of the entanglement entropy exhibit distinct features after quantum quenches in the anisotropic next-nearest neighbor Ising (ANNNI) chain upon tuning the quench parameter across an underlying quantum phase transition. Specifically, we find that the derivatives of local observables with respect to the quench parameter develop prominent dips/peaks in the vicinity of the quantum phase transition. We determine the quantum real-time evolution by means of the infinite-Time Evolved Block Decimation (iTEBD), which provides numerically exact results for the transient to intermediatetime dynamics in the thermodynamic limit [10][11][12] .
In order to access the long-time (asymptotic) properties of the considered quantum-chaotic system after the expected thermalization, we employ a numerical linked cluster expansion (NLCE) for thermal equilibrium states 13 . We again find distinct signatures of the quantum phase transition in derivatives of the correlation functions. Also, the effective temperature exhibits a marked minimum as function of the quench parameter in the close vicinity of the quantum phase transition. Since the considered one-dimensional system does not support singular behavior after equilibration, upon assuming that eigenstate thermalization occurs [14][15][16][17] , these prominent features are not associated with nonanalytic properties (in contrast to the integrable systems studied in Refs. 6,7 ), but nevertheless represent distinct signatures of quantum phase transitions. Finally, we discuss similar phenomena for quantum phase transitions involving topologically different quantum states. We also discuss how our findings can be tested in current experiments with Rydberg atoms.
The presentation is organized as follows. In Sec. II, we introduce the Hamiltonian of the ANNNI chain, and introduce the protocol used to probe the ferromagnetic to paramagnetic quantum phase transition via dynamics following quantum quenches. The results obtained for dynamics after the quenches are presented in Sec. III, while the results after thermalization are presented in Sec. IV. Combining results from the dynamics and thermalization, in Sec. V we report the estimated phase diagram for the ferromagnetic to paramagnetic quantum phase transition for a wide range of parameters of the ANNNI chain. In Sec. VI, we discuss the feasibility of testing our results experimentally, while in Sec. VII we discuss the applicability of our protocol to detect topological quantum phase transitions. In Sec. VIII, we summarize our results and discuss their implications.

II. ANNNI HAMILTONIAN AND QUENCH PROTOCOL
The ANNNI chain is a very well studied spin model (see, e.g., Ref. 18 ). Its Hamiltonian in a chain with L sites can be written aŝ When mapped onto a fermionic Hamiltonian using Jordan-Wigner transformation 19 , the next-nearestneighbor term (with strength κ) maps onto a fourfermion interaction. At T = 0, this model has a rich (and still partly controversial) phase diagram in the κ-Γ plane. The quantum phase transition line from the ferromagnetic to the paramagnetic phase, which occurs as the antiferromagnetic next-nearest neighbor coupling κ > 0 crosses a critical value for a fixed |Γ| < 1, is a second order phase transition (see, e.g., Ref. 20 ). In the quadrant (κ > 0, Γ > 0), this line is well described using second order perturbation theory, with the critical parameters satisfying (see Fig. 11) 20 : To probe this ferromagnetic to paramagnetic quantum phase transition at a fixed value of Γ, we generate a family of HamiltoniansĤ(κ). We then generate a family of nonequilibrium states via quenches withĤ(κ). The protocol (straightforward to generalize to other models) consists of following steps (see Fig. 1): (i) The initial state is fixed to be the ground state |ψ(κ I ) ofĤ(κ I ), where κ I is deep in the ferromagnetic phase.
(iii) We compute expectation values of observables O(t, κ) = ψ(t, κ)|Ô|ψ(t, κ) . (iv) For a fixed value of t, we study how O(t, κ) changes with κ, focusing on the behavior in the vicinity of κ c , where κ c is the critical value of κ for the transition given the selected value of Γ.

III. QUANTUM DYNAMICS
We study the time evolution of observables after quantum quenches in an infinite ANNNI chain using iTEBD (see Appendix IX) [10][11][12] . Following the protocol introduced in Sec. II (see Fig. 1), we fix Γ (we take Γ = 0.2) and then fix κ I so that the initial state is a ground state of the ANNNI chain deep in the ferromagnetic phase (we take κ I = 0). In our iTEBD calculations, we introduce a very small (∼ 10 −6 ) longitudinal field to pin one of the two degenerate maximally polarized ground states. The critical value of κ for the ferromagnetic to paramagnetic quantum phase transition for Γ = 0.2 is κ c ≈ 0.41.
We first focus on the dynamics of two local observables, the nearest and next nearest neighbor longitudinal correlators In for six values of κ after the quench. The dynamics of both longitudinal correlators is qualitatively similar for the values of κ shown. Their decrease with time speeds up as κ increases about κ c . How the closeness to κ c affects the dynamics is better seen by plotting the correlations for fixed times t after the quench as functions of κ [step (iv) in the protocol introduced in Sec. II]. This is done in Fig. 3, where we show results for C x 1 [ Fig. 3(a)] and C x 2 [ Fig. 3(b)]. At all times reported, C x 1 and C x 2 decrease rapidly with increasing the value of κ for κ κ c . In addition, with increasing time, the decrease in the correlators becomes more prominent when κ ≈ κ c . This is apparent in the insets, where we show the derivative of the correlators. They develop sharper dips close to κ c as the evolution time increases.
In Refs. 6,7 it was proved that following the same protocol discussed here but for noninteracting models (or models mappable to them) results in nonanalytic behavior of local observables at the quantum phase transition in the limit t → ∞ (after having taken the thermodynamic limit first). While this is not the case in the quenches in generic models studied here (see Sec. IV), the prominent features seen in Fig. 3 at finite times are promising for an experimental determination of κ c .
We also studied the dynamics of the half chain entanglement entropy S 1/2 = −Tr[ρ 1/2 ln ρ 1/2 ], where ρ 1/2 is the density matrix of the half chain (obtained by tracing out the other half). This is a nonlocal observable that is expected to increase linearly with time in quantumchaotic systems 21 . In Fig. 4(a), we plot the time evolution of S 1/2 for six values of κ. As for the local operators in Fig. 2(a), the change of S 1/2 with time speeds up as κ increases about κ c . Figure 4(b) shows S 1/2 at fixed times t after the quench plotted vs κ, and the inset in Fig. 4(b) shows the derivative with respect to κ of the results in the main panel. Like the local operators in Fig. 3, the behavior of the half chain entanglement entropy carries a marker of the quantum phase transition.

A. Changing κI
In Fig. 5, we show the derivative of C x 1 with respect to κ at a fixed time after the quench, plotted as a function of κ, for different values of κ I in the initial ground state. We recall that as κ I departs from κ c , for κ I < κ c , the ground state of the system is deeper in the ferromagnetic phase. The results in Fig. 5 show that starting deeper in the ferromagnetic phase results in a slightly shallower dip in dC x 1 /dκ, while its position remains unchanged. This might be expected as the departure of κ I from κ c increases the magnitude of the quench, and hence increases the final energy density, thereby blunting the signature of the quantum phase transition. This is consistent with the results in Sec. IV C, where we discuss the effect that the increase in the magnitude of the quench has in observables after thermalization.

IV. RESULTS AFTER THERMALIZATION
Because of the linear growth of the entanglement entropy seen in Fig. 4, the iTEBD technique only allows one to study dynamics at short and intermediate times. To explore the fate of observables after thermalization, we use a numerical linked cluster expansion (NLCE) 13 . We broaden the class of initial states to explore how initial nonzero temperatures modify the behavior of observables after thermalization.
Here we consider more general quenches within the ANNNI Hamiltonian involving initial statesρ I that are thermal equilibrium states of the initial Hamiltonian H(κ I ). For an initial temperature T I ,ρ I has the form When T I = 0,ρ I is the ground state ofĤ(κ I ). As in the previous section, we quench κ I → κ, while Γ is kept unchanged (Γ = 0.2). In Secs. IV A and IV B we fix κ I = 0. In Sec. IV C, we explore what changes when κ I is varied within the ferromagnetic phase (κ I < κ c ≈ 0.41). Since the energy after the quench is the only conserved quantity, at sufficiently long times in the thermodynamic limit, observables are expected to be described by a Gibbs ensemble 17ρ with a temperature T (κ) > 0 (which is nonzero even when T I = 0) determined by the energy E(κ) set by the initial stateρ I , as dictated by: We use the numerical linked cluster expansion (NLCE) technique introduced in Refs. 13 to study the thermal expectation values of observables in the thermodynamic limit (see Appendix IX for details). All the NLCE results for the ANNNI chain are obtained using 15 orders of the maximally connected cluster expansion introduced in Refs. 22 . To gauge how well the series has converged, we estimate the convergence error for an observable by computing the relative difference between the last two orders (14 and 15) of the NLCE 22 . We only report results whose convergence error for the energy is less than 10 −5 . T (κ) is obtained by numerically matching the energies in the left and right side of Eq. (6). Both energies are evaluated using NLCE to 15 orders, and T (κ) is computed by enforcing that their relative difference is less than 10 −11 (see Ref. 22 ). For observables other than the energy, we only report results whose convergence errors are less than 5 × 10 −5 (except for the entropy, for which we set the cut off to be 7 × 10 −5 ). Those errors are small enough to be unimportant for the discussions that follow.

A. Observables
As mentioned before, in the thermodynamic limit at sufficiently long times after the quench, thermalization is expected to occur in the nonintegrable systems considered here 17 . Next we study the expected thermal equilibrium results that observablesÔ reach after equilibration following the quench.
In the space of all possible thermal equilibrium ensembles parameterized by the coordinates (T, κ), the initial stateρ I sets a trajectory T (κ) determined by Eq. (6). One can then write Since O(T, κ) is an analytic function whenever T > 0, and since dT /dκ is expected to be a smooth function of κ (we discuss this in Sec. IV B), then dO/dκ must be −10 (3), evaluated in thermal equilibrium using NLCE following quenches κI = 0 → κ, with TI = 0, 0.1, 0.5, and 1.0. We also show C x 1 and C x 2 in the ground state ofĤ(κ) (dotted lines) computed using iTEBD. The main panels in (a) and (b) show dC1/dκ and dC2/dκ, respectively, while the corresponding insets show C x 1 and C x 2 . The vertical dashed lines mark the critical κc ≈ 0.41. a smooth function of κ after equilibration following the quench. Still, for observables that are indicators of the quantum phase transition in nonintegrable systems (e.g., order parameters and related observables), (∂O/∂T ) κ and (∂O/∂κ) T can be large if T is low when κ is close to κ c (we show the latter to be the case for our quenches in Sec. IV B). This means that, even in thermal equilibrium, it is possible to have prominent (but smooth) features in dO/dκ as observed at intermediate times in the previous section. In integrable systems, in which all possible states after equilibration are described by generalized Gibbs ensembles that are parameterized by extensive numbers of quantities 23,24 , nonanalytic behavior is possible and has in fact been observed in Refs. 6,7 .
In Fig. 6, we show the thermal equilibrium results obtained for the nearest C x 1 and next nearest C x 2 neighbor longitudinal spin correlations per site [see Eq. (3)] as functions of κ after the quench, as well as their expectation values in the ground state ofĤ(κ) computed with iTEBD. The main panels show dC 1(2) /dκ, while the insets show C 1(2) (κ), for various initial temperatures T I and in the ground state (dotted lines, computed with iTEBD). In the ground state, C x 1 and C x 2 are nearly one in the ferromagnetic phase and exhibit a rapid decrease when crossing the quantum phase transition (prominent minima can be seen in dC 1(2) /dκ at κ c ), i.e., they serve as indicators of the quantum phase transition. (They also serve as indicators of the ferromagnetic to paramagnetic quantum phase transition in the integrable transverse field Ising model, see Appendix X.) This zerotemperature behavior is the precursor of the behavior of C x 1 and C x 2 observed in the insets for low initial T I , which, in turn, produces the prominent minima in dC 1(2) /dκ near κ c observed in the main panels. Figure 6 shows that the position of the minima drift away from κ c , and they become shallower, with increasing T I . Note that the results for T I = 0 and T I = 0.1 overlap in the plots.
Qualitatively similar results were obtained for other local observables, such as the transverse magnetization per site m z = i σ z i /L, shown in Fig. 7(a), and for the (von-Neumann) entropy per site s = −tr(ρ GE lnρ GE )/L of the thermal stateρ GE (κ), shown in Fig. 7(b). In the ground state, m z increases rapidly when transitioning from the ferromagnetic to the paramagnetic phase, as shown in Fig. 7(a) (dotted lines, computed with iTEBD). Hence, m z serves as an indicator of the quantum phase transition, and its behavior at zero temperature is the reason there are prominent maxima in dm z /dκ near κ c for quenches at low T I . (See Appendix X for ground-state results of m z across the ferromagnetic to paramagnetic quantum phase transition in the integrable transverse field Ising model.) The entropy, on the other hand, is strictly zero at zero temperature, i.e., it does not change at the phase transition [the entanglement entropy does change, as shown in Fig. 4(b)]. However, as we show in Sec. IV B, when T I is low, the temperature after quench increases rapidly when κ crosses κ c and this produces the rapid increase of s seen in Fig. 7(b).

B. Temperature
Let us now show that dT /dκ is a smooth function of κ. The ANNNI Hamiltonian can be written asĤ(κ) = H 0 + κV , so that keeping the initial stateρ I fixed and changing κ after the quench results in E(κ) being a linear function of κ As in the previous section for dO/dκ, for the energy one can write where (∂E/∂T ) κ = C κ (T ) is the specific heat. Combining Eqs. (8) and (9), we have that All functions in the r.h.s. of Eq. (10) are smooth, and C κ [T (κ)] > 0, because T (κ) > 0 after the quench. This shows that T (κ) is also a smooth function. Next, we use numerical calculations to explore whether quenches κ I → κ spanning across κ c produce temperatures T (κ) with signatures of the quantum phase transition, as shown to be the case in Sec. IV A for local observables. Figure 8 shows T (κ) for quenches with κ I = 0 → κ for various initial temperatures T I , including the ground state ofĤ(κ I ). For very low initial temperatures T I 0.1, the temperatures T (κ) after the quench are essentially indistinguishable from those for T I = 0. This explains why all the results reported in Sec. IV A are indistinguishable for T I = 0 and T I = 0.1. For those very low T I , the temperatures T (κ) exhibit a low-temperature minimum in the vicinity of κ c (at κ m ≈ 0.39, for which T (κ m ) ≈ 0.06). At T I = 0.5, a temperature at which T (κ) after the quench departs from the T I = 0 result, a minimum in T (κ) still remains visible close to κ c . The locus of minima in T (κ), shown as a dotted line for a large number of T I , makes apparent that the minima remain close to κ c as long as T I remains low (T I 1.0). At higher initial temperatures, the minima depart from κ c indicating that the information about κ c is washed out.
Overall it is remarkable that, due to the presence of the phase transition (and the corresponding closing of the gap above the ground state), when quenching to the where we defined the intensive counterparts of the extensive quantities in Eq. (10) as a = A/L and e = E/L. The main panel in Fig. 9 shows (∂e/∂κ) T vs κ at different temperatures [inset Fig. 9(a) shows e vs κ at the same temperatures]. For T = 0, we also show iTEBD results (the NLCE results do not converge close to κ = κ c ). Notice that, in the region in which the NLCE results converge to the precision mentioned in the introduction of this section, they are indistinguishable from the iTEBD ones. The iTEBD results for (∂e/∂κ) T =0 exhibit a rapid decrease about κ c [resulting in a singularity in (∂ 2 e/∂κ 2 ) T =0 at κ c , as shown in inset Fig. 9(b)], reflecting the nonanalytic behavior of the energy at the (second order) quantum phase transition. That rapid decrease leaves its signature in the low-temperature behavior of (∂e/∂κ) T >0 , and this is what makes possible for Eq. (11) to be satisfied close to κ c for low initial temperatures.
In Fig. 9, a = σ x i σ x i+2 ρ I /L is shown as a horizontal line for T I 0.1. For those very low initial temperatures, a is very close to 1 (a ≈ 0.99) since κ I = 0 is deep in the ferromagnetic phase, and Fig. 9 shows that the condition (∂e/∂κ) T = a is satisfied at κ = 0.39 for T = 0.05 and at κ = 0.37 for T = 0.1. Those two temperatures approximately bound the range of effective temperatures after the quench for κ close to κ c when T I 0.1, see Fig. 8. This explains why the minimum in T (κ) vs κ occurs very close to κ c for T I 0.1. Increasing the initial temperature beyond T I = 0.1 increases T but also reduces the value of a. This results in the minimum remaining close (and actually slightly approaching) κ c in Fig. 8 when T I departs from 0.1 but still remains low (T I 1.0). Since the slope of (∂e/∂κ) T at the crossing point near κ c is negative, it follows from Eq. (10) that the extremum in T (κ) near κ c is a minimum.

C. Changing κI
Motivated by the results discussed in Sec. III A, we explore next what happens to the thermal equilibrium results after equilibration when one changes κ I within the ferromagnetic regime, keeping T I = 0 fixed. In Fig. 10(a), we show T (κ) vs κ for κ I = −0.2, 0, and 0.2. As expected from the fact that the initial state remains a nearly perfect ferromagnet, the minima in T (κ) close κ c are robust to the choice of initial κ I . However the minimum value of T (κ) attained decreases as κ I approaches κ c . As a result, the signature of the presence of a quantum critical point in observables after thermalization becomes sharper as κ I → κ c . This is apparent in Fig. 10(b) in which we plot dC x 1 /dκ. Note that in Fig. 10(a) there is a singularity in T (κ) at κ = 0.2 for κ I = 0.2, as well as at κ = 0 for κ I = 0. These are trivial consequences of performing no quench, which means that the system remains in the ground state. The fact that |dT /dκ| → ∞ at those points follows from Eq. (10) due to specific heat C κ (T → 0) → 0 in the denominator. These singularities have no consequence in the expectation values of observables.

V. PHASE DIAGRAM
Here we combine results obtained for C x 1 at intermediate times after the quench (from iTEBD calculations), and after equilibration (from NLCE calculations), to identify, in the (κ, Γ) plane, the phase boundary separating the ferromagnetic and paramagnetic phases in the ground state. We estimate κ c by carrying out quenches κ I = 0 → κ for different values of Γ (Γ is not changed during the quench). Qualitatively similar results were obtained for other local observables such as C x 2 and m z and are not reported here.
In the main panel of Fig. 11, we show κ c extracted from the extrema of dC x 1 /dκ obtained using iTEBD results at t = 25 after quenches starting from the ground state, and NLCE thermal equilibrium results after quenches starting from the ground state (T I = 0) and from an initial temperature T I = 0.3. As Γ increases, the NLCE convergence errors are higher for quenches starting from the ground state. This occurs because the critical point gets closer to κ I = 0 and the effective temperature after the quench becomes too small (see Fig. 10 and related discussion). This is the reason no NLCE points are reported for quenches with Γ ≥ 0. critical point. This also affects the NLCE convergence, resulting in no NLCE data points for Γ 0.2. The results in Fig. 11 show that both the intermediate-time and (expected) long-time extrema follow very closely the phase boundary calculated using iTEBD for the ground state [locating the singularity in (∂ 2 e/∂κ 2 ) T =0 ], which is well described by the second order perturbation theory results.

VI. EXPERIMENTAL TESTS
It is a central aspect of this work that the reported signatures of the quantum phase transitions in the ANNNI model are accessible in state-of-the-art quantum simulator platforms with Rydberg atoms. The ANNNI Hamiltonian (1) can be straightforwardly realized using using Rydberg dressing in ultracold atoms in optical lattices 25,26 . Rydberg-dressed atoms exhibit a softcore interaction potential J i,j = J 0 /[1 + (R ij /R c ) 6 ], which is approximately constant below a threshold distance R c between two atoms and decays quickly beyond the threshold R c (in a R −6 ij fashion as a function of distance R ij = a|i − j|, where a is the lattice spacing between the involved spins) 25 . Realizing approximately the ANNNI model with such a soft-core interaction potential requires to choose the tunable threshold R c such that (R i,i+3 /R c ) 6 ( In such a regime only nearest and next-nearest neighbor couplings have to be taken into account, while further distant ones can be neglected. The relative strength of nearest and next-nearest neighbor interactions, quantified by κ = J i,i+2 /J i,i+1 = [1 + (a/R c ) 6 ]/[1 + (2a/R c ) 6 ] in Eq. (1), can also be varied by tuning R c relative to the lattice spacing a with the only limitation that κ < 1. As the targeted quantum critical point κ c ≈ 0.41 < 1, the reported signatures therefore lie within the tunability of the couplings. Let us note that the interaction in the experiment would be of antiferromagnetic nature and not directly of the type required in Eq. (1). However, by performing a rotation σ x l → −σ x l on every other lattice site, e.g., even ones, the Hamiltonian in Eq. (1) maps onto a purely antiferromagnetic spin model and therefore to the one which can be realized experimentally. Furthermore, transverse fields can be straightforwardly generated, implying that the full Hamiltonian can be modeled with high accuracy.
It remains to clarify whether also the dynamics of this system can be accessed in the desired regimes, which we now answer in the affirmative. First, Rydberg-dressed atom systems with a large number of spins (L ≈ 200) were already created in Ref. 25 . The trapping potential for the ultracold quantum gas only has a minor impact when considering Rydberg dressing, it affects the preparation of the initial condition by limiting the maximal number of spins which can be controllably initialized 25 . Specifically, the fully polarized initial condition we are considering in our work can be prepared with high fidelity as demonstrated in Ref. 26 . Hence, the main point that remains to be addressed is the coherence time, i.e., whether it is possible to identify the proposed signatures before decoherence sets in. In a recent experiment with Rydbergdressed atoms time scales Jt 10 were achieved, where J denotes the strength of the nearest-neighbor couplings. Consequently, the time scales discussed in Sec. III are in the experimentally accessible regime. We note that also the desired spin-spin correlation functions in Eq. (3) can be measured in the aforementioned experimental systems 26 .

VII. TOPOLOGICAL TRANSITIONS
A final question we address next is how generally one can use the previously introduced protocol to locate quantum phase transitions in one-dimensional models. Given the results obtained and insights gained within the ANNNI chain (notice that in Fig. 11 we report results for an entire phase boundary), we expect this protocol to be widely applicable to one-dimensional models with traditional quantum phase transitions. A different question is whether such signatures in local quantities can be used to locate topological quantum phase transitions, as shown for noninteracting models in Ref. 7 (non-local quantities can, of course, retain such information in the noninter- acting case -see, e.g., 27,28 ). In what follows, we report results from a preliminary exploration of dynamics after quantum quenches about topological transitions in two quantum-chaotic models.
First, we explore the quantum phase transition from the Néel to the symmetry protected topological "Haldane" phase in the spin-1 anisotropic (XXZ) Heisenberg chain model. The Hamiltonian for this model readŝ whereŜ x,y,z i denote the x, y and z components of the spin-1 operator at site i. Four different phases occur in this model when one changes the anisotropy parameter ∆ (see, e.g., Refs. 29,30 and references therein). Here we focus on the transition that occurs upon decreasing ∆ from ∆ > 1, a limit in whichĤ XXZ reduces to the spin-1 Ising antiferromagnet. With decreasing ∆, the ground state of H XXZ undergoes a quantum phase transition from the antiferromagnet to the Haldane phase at ∆ c ≈ 1.183. The Haldane phase is a topological phase, protected by any one of the following three global symmetries: D 2 spin rotation, time-reversal, and bond centered inversion 31 . This transition is of second order, and belongs to the 2D Ising universality class 32,33 .
In Fig. 12, we show ground state results for dC x 1 /d∆, and for dC z 1 /d∆, where plotted as functions of ∆. As for the local observables shown in Fig. 6 for the ANNNI model, dC x 1 /d∆ in Fig. 12(a) [dC z 1 /d∆ in Fig. 12(b)] exhibits a sharp maximum (minimum) at the transition point. We expect this maximum (minimum) to be the precursor of a peak (dip) close to ∆ c after the quantum dynamics generated following the protocol introduced in Sec. II. To test this, we take as initial state the ground state at large ∆ I = 2 and quench ∆ across the neighborhood of ∆ c . Due to the high computational cost of the iTEBD calculations for the spin-1 anisotropic Heisenberg chain, we are only able to study dynamics at short times (t ≤ 7) after the quench. Still, for these short times, Fig. 12(a) [ Fig. 12 shows that a peak (dip) appears to develop in dC x 1 /d∆ (dC z 1 /d∆) about a ∆ * greater than, but close to, the transition point ∆ c . As t increases those peaks sharpen and move toward ∆ c . This suggest that our protocol can be used to locate this phase transition.
In the spin-1 anisotropic (XXZ) Heisenberg chain model, like in the ANNNI model, to locate the phase transition using our protocol we rely on the rapid change of local correlations close to the transition point. Next, we study a model in which at the transition point in equilibrium, due to a symmetry, there is a vanishing change in local correlations. The question then is whether this can also be detected in the quantum dynamics and used to locate the transition point.
The model is the bond-alternating Heisenberg model 34 where σ i are the Pauli matrices (periodic boundary condition implied). This model exhibits a topological transition between two dimerized phases at η c = 1, which can be located using a nonlocal string order parameter. Because of the invariance (up to a rescaling) of Hamiltonian (15) under η → 1/η, one can see that in thermal equilibrium local correlations are symmetric about η c = 1. This means that, so long as the correlations depend on η, they must exhibit a maximum or a minimum at η c = 1. In Fig. 13 we show that this is indeed the case for C x 1 and C x 2 , defined in Eq. (3), in the ground state. At η c = 1, C x 1 exhibits a minimum in Fig. 13(a), and C x 2 exhibits a maximum in Fig. 13(b). Next, we explore the fate of those extrema in the quantum dynamics. We quench the parameter η following our protocol in Sec. II, namely, taking the initial state to be the ground  Figure 13 shows that extrema occur at η * close to η c , both in the short-time dynamics (studied using iTEBD, and shown in the main panels) and after thermalization (studied using NLCE, and shown in the insets). We note that the position η * of the minimum in C x 1 (maximum in C x 2 ) relative to η c depends on whether the initial state chosen is the ground state for η I greater or smaller than η c . The minimum (maximum) develops at η * < η c for η I < η c , as seen in Fig. 13(a) [Fig. 13(b)], and would develop at η * > η c for η I > η c . This is also a result of the invariance (up to a rescaling) of Hamiltonian (15) under η → 1/η. Hence, our protocol can also be used in this model for which the transition is located directly using the local observables, as opposed to the earlier models for which it was located using derivatives of the local observables.

VIII. SUMMARY AND DISCUSSION
In summary, we have shown that local observables can be used to locate the ferromagnetic to paramagnetic quantum phase transition in the ANNNI chain (a nonintegrable model) both at intermediate times after a quench and at long times after thermalization. The initial states for our quenches were chosen to be ground states of the ANNNI chain deep in the ferromagnetic phase. We explored the effect that changing the magnitude of the quench and starting from initial finite temperature states has in many of our conclusions, and showed that our conclusions are robust against those changes. We also discussed potential experimental tests, as well as the applicability of our protocol to detect topological quantum phase transitions.
More generally, the fact that intermediate-time dynamics, following quenches whose initial states are ground states far from a quantum phase transition, provide a way to locate the quantum phase transition is promising for experiments with ultracold quantum gases 3,4 and ions 35,36 . In those experiments, it is usually straightforward to prepare ground states far away from quantum phase transitions but it is much more challenging to prepare them close to the transitions. The latter is needed to locate the quantum critical point via traditional measurements of the system in equilibrium. Also, not needing to wait long times to observe signatures of the quantum phase transition in the dynamics after the quench is important because, due to heating and other undesirable effects, keeping the dynamics coherent in the experiments becomes increasingly challenging as the evolution time increases. support by the Deutsche Forschungsgemeinschaft via the Gottfried Wilhelm Leibniz Prize program. This research was also supported in part by the International Centre for Theoretical Sciences (ICTS) during a visit for the program -Thermalization, Many body localization and Hydrodynamics (Code: ICTS/hydrodynamics2019/11).

A. Infinite Time-Evolving Block Decimation
In this section, we briefly outline details about the infinite-time evolving block decimation (iTEBD) method. The iTEBD algorithm is based on the infinite matrix-product state (iMPS) representation, which can efficiently represent many-body wave functions with the accuracy controlled by the bond dimension χ (the error decreases rapidly with increasing χ). A general quantum state |Ψ on a chain with L sites can be written in the following MPS form 10,11 : where, A[n] sn is a χ n−1 ×χ n dimensional matrix and |s n with s n = 1, . . . , d is a basis of local states at site n. Doing repeated Schmidt decomposition on the state |Ψ , one can get the form for the coefficients c s1...s L : where Γ s are rank-3 tensors, and Λ s are positive, real, square-diagonal matrices. After doing the tensor contractions, the structure obtained can be readily identified with a Matrix Product State as in equation (16).
The size of the tensors χ i required to represent a state can be shown to be related to the von Neumann entropy S i of the partition 1 . . . i : i + 1 . . . L, as S i ≤ 2 ln χ i . If the entropy is area-law (as is the case for ground states of one-dimensional gapped systems), χ i remains finite in the thermodynamic limit.
Using the iTEBD algorithm, one can evaluate the time evolution of a quantum state: and use the imaginary time evolutionÛ (τ ) = exp(−Ĥτ ) to find the ground state of the HamiltonianĤ. Using the Trotter-Suzuki decomposition to the first order, one can write whereÂ andB are operators, and δ is a small parameter. and decompose it as a sum The terms within one partition act on different sites and thus commute with each other: One can approximate the time evolution operator for a very small time slice δt 1, to the first order, using (19), as: To determine the suitable δ, one can successively make it smaller to achieve convergence. We used the bond-link dimension χ = 4000 to ensure convergence for the longest real-time dynamics results and the time steps used is δ = 0.01. The time evolution in equation (18) is obtained by applying the operators e −iĤ odd δ and e −iĤevenδ iteratively to the initial state |ψ(0) , which has been previously decomposed in the form of an MPS. After the application of each operator at sites i and i + 1 the decomposition (17) (17), where Γ [i] and λ [i] are independent of i. Thus, given that the time evolution is generated by two-site operators, only the tensors Γ A , Γ B , λ A , and λ B have to be updated, where Γ A = Γ [2i] , Γ B = Γ [2i+1] , λ A = λ [2i] , and λ B = λ [2i+1] .
In our case, in which we also have a next-nearest neighbor interaction, one can group the sites (merge two neighboring site to one) and proceed with the same algorithm where the local Hamiltonian is now 16×16 instead of 4×4.

B. Numerical linked cluster expansion
For lattice models in the thermodynamic limit (L → ∞), NLCE allows one calculate the expectation value of extensive observablesÔ per site, O = Ô /L, as a sum over contributions from all connected clusters c that can be embedded on the lattice: where W O (c) is the weight of cluster c, and M (c) is the number of ways per site in which one can embed c on the lattice. W O (c) is computed for each cluster c using the inclusion exclusion principle: where Ô c is the expectation value ofÔ in the cluster c, and the sum runs over all connected sub-clusters of c.
For the smallest cluster c 0 , W O (c 0 ) = Ô c0 . For each cluster, Ô c = Tr[ρ cÔ ], whereρ c is the relevant density matrix in the cluster. For the initial stateρ c is of the form Eq. (4), and for the thermal state used to describe observables after equilibrationρ c is of the form Eq. (5), with their respective Hamiltonians restricted to the cluster c. Ô c is calculated numerically using full exact diagonalization.
We use the maximally connected expansion introduced in Ref. 22 , in which each cluster c contains all possible bonds between the sites as per the specific Hamiltonian considered. The order of the NLCE is then the number of lattice sites of the largest cluster c considered in the sum (23). The series is convergent when errors in consecutive orders vanish exponentially fast with increasing order.
For the thermal equilibrium results in the bondalternating Heisenberg model in Sec. VII [ Fig. 13], we calculate Ô c separately for the bond-alternating Hamiltonian and its reflected configuration in c and average them. This extra step is necessary to restore the translational invariance assumed to build the NLCE used. For this model, in order to calculate the temperatures of the thermal equilibrium ensembles used to describe observables after thermalization, we use energies after the quench that are obtained using iTEBD. With those energies, the temperatures are obtained using a 16-order NLCE calculation. The convergence errors in the calculation of the energy are smaller than 5 × 10 −4 for all parameters considered (they are much smaller than 5 × 10 −4 for most parameters considered).

X. TRANSVERSE-FIELD ISING CHAIN
The transverse field Ising chain (TFIM) is probably the most studied exactly solvable (integrable) model in the context of quantum phase transitions 1,18 . Its Hamiltonian readsĤ It is the noninteracting limit (κ = 0) of our ANNNI Hamiltonian [Eq. (1)]. In Fig. 14, we report ground state results for C 1(2) and m z [ Fig. 14(a)], and their derivatives [ Fig. 14(b)], across the ferromagnetic to paramagnetic phase transition, which occurs in this model at Γ = 1.