Determination of dynamical quantum phase transitions in strongly correlated many-body systems using Loschmidt cumulants

Dynamical phase transitions extend the notion of criticality to non-stationary settings and are characterized by sudden changes in the macroscopic properties of time-evolving quantum systems. Investigations of dynamical phase transitions combine aspects of symmetry, topology, and non-equilibrium physics, however, progress has been hindered by the notorious difficulties of predicting the time evolution of large, interacting quantum systems. Here, we tackle this outstanding problem by determining the critical times of interacting many-body systems after a quench using Loschmidt cumulants. Specifically, we investigate dynamical topological phase transitions in the interacting Kitaev chain and in the spin-1 Heisenberg chain. To this end, we map out the thermodynamic lines of complex times, where the Loschmidt amplitude vanishes, and identify the intersections with the imaginary axis, which yield the real critical times after a quench. For the Kitaev chain, we can accurately predict how the critical behavior is affected by strong interactions, which gradually shift the time at which a dynamical phase transition occurs. We also discuss the experimental perspectives of predicting the first critical time of a quantum many-body system by measuring the energy fluctuations in the initial state, and we describe the prospects of implementing our method on a near-term quantum computer with a limited number of qubits. Our work demonstrates that Loschmidt cumulants are a powerful tool to unravel the far-from-equilibrium dynamics of strongly correlated many-body systems, and our approach can immediately be applied in higher dimensions.


I. INTRODUCTION
Whether or not quantum many-body systems out of equilibrium can be understood in terms of well-defined phases of matter is a central question in condensed matter physics. The lack of universal principles, such as those governing equilibrium systems [1,2], makes the problem exceptionally hard. Still, the concepts of criticality and far-from-equilibrium dynamics have recently been elegantly unified through the discovery of dynamical phase transitions in which a time-evolving quantum many-body system displays sudden changes of its macroscopic properties [3][4][5][6][7][8][9][10][11]. In equilibrium physics, phase transitions are reflected by singularities in the free energy, and dynamical phase transitions are similarly given by critical times, where a non-equilibrium analogue of the free energy becomes non-analytic. Specifically, the role of the partition function is played by the return, or Loschmidt, amplitude of the many-body system after a quench, and its logarithm yields the corresponding free energy.
A typical setup for observing dynamical quantum phase transitions is depicted in Fig. 1a: a onedimensional chain of interacting quantum spins is initialized in a ground state characterized by one type of order (or the lack of it) and subsequently made to evolve according to a Hamiltonian whose ground state possesses a different order. Experimentally, dynamical phase transitions have been observed in strings of * Email: teemu.ojanen@tuni.fi trapped ions [12,13], optical lattices [14], and several other systems that offer a high degree of control [15][16][17][18][19]. The Loschmidt amplitude is the overlap between the initial state of the system and the state of the system at a later time. Moreover, similarly to equilibrium systems, dynamical phase transitions may occur at critical times, where the Loschmidt amplitude vanishes, and the dynamical free energy becomes non-analytic in the thermodynamic limit. As illustrated in Fig. 1b, these non-analytic signatures may appear as cusps in the dynamical free energy, however, strictly speaking, they only occur for infinitely large systems. For finite-size systems, they are typically smeared out, and often for spin chains at least L 50 − 100 spins are required in order to identify and determine the critical times of a dynamical phase transition. Since the Hilbert space dimension grows exponentially with the chain length, the outstanding bottleneck for theoretical investigations of dynamical phase transitions is the need to predict the far-from-equilibrium dynamics of large quantum systems. Numerically, the task is computationally costly, or even intractable, and generally it requires advanced system-specific techniques that do not easily generalize to other systems or spatial dimensions [4,5,10,[20][21][22][23][24][25][26][27][28][29][30][31][32][33]. For this reason, little is still known about dynamical phase transitions and the general applicability of concepts like universality and scaling. In fact, our current understanding comes to a large extent from a few exactly solvable models [3-11, 26, 34-43]. Important questions concern the relationship between critical times and dynamical changes in local observables or the entanglement spectrum (or other dynamical measures), which often exhibit similar but not strictly re- Dynamical phase transitions. a, A sudden quench of the system parameters causes a dynamical phase transition in a quantum spin chain with L sites. b, In the thermodynamic limit, such phase transitions give rise to singularities in the rate function at the critical times, tc,1, tc,2 . . ., see Eqs. (1) and (2) for definitions, however, in finite-size systems they are smeared out. c, The singularities in the rate function are associated with the zeros (circles) of the Loschmidt amplitude in the complex-time plane. In the thermodynamic limit, they form continuous lines, and the real critical times are given by the crossing points with the imaginary axis. We determine the zeros of the Loschmidt amplitude from the Loschmidt cumulants evaluated at the basepoint τ . lated time scales. However, case-by-case investigations have revealed that dynamical phase transitions are often accompanied by interesting dynamics with comparable time scales, and one could view them as indicators of nontrivial dynamics in other many-body properties.
Here we pave the way for systematic investigations of dynamical phase transitions in correlated systems using Loschmidt cumulants, which allow us to accurately predict the critical times of a quantum many-body system using remarkably small system sizes, on the order of L 10 − 20. Using modest computational power, we determine the critical times of the interacting Kitaev chain and the spin-1 Heisenberg chain after a quench and find for instance that a dynamical phase transition in the Kitaev chain gets suppressed with increasing interaction strength. The Loschmidt cumulants allow us to determine the complex zeros of the Loschmidt amplitude as illustrated in Fig. 1c. We can thereby map out the thermodynamic lines of zeros and identify the crossing points with the imaginary axis, corresponding to the real critical times, where a dynamical phase transition occurs. This approach makes it possible to predict the critical dynamics of a wide range of strongly interacting quantum many-body systems and is applicable also in higher dimensions. In two dimensions, the zeros can make up lines or surfaces in the complex plane, and our method can be used to determine all of these zeros as well as their density. Moreover, as we will show, our method provides exciting perspectives for future experiments on dynamical phase transitions. Specifically, our method makes it possible to predict the first critical time of a quantum many-body system after a quench by measuring the fluctuations of the energy in the initial state. In addition, due to the favorable scaling properties of our method, it is conceivable that it can be implemented on a near-term quantum computer with a limited number of qubits.
We now proceed as follows. In Sec. II, we develop our method for determining the zeros of the Loschmidt echo and their crossing points with the imaginary axis in the thermodynamic limit, which yield the critical times of a quantum many-body system after a quench. In Sec. III, we consider dynamical phase transitions in the Kitaev chain after a quench, and we show how we can determine the critical times from remarkably small chain lengths even with strong interactions. Section IV is devoted to the spin-1 Heisenberg chain and includes several quenches, for instance, from the Haldane phase to the Néel phase. In Sec. V, we discuss the experimental perspectives for future realizations of our method. Finally, in Sec. VI, we state our conclusions and provide an outlook on possible avenues for further developments.

II. FROM LOSCHMIDT CUMULANTS TO LOSCHMIDT ZEROS
The fundamental object that describes dynamical phase transitions is the Loschmidt amplitude, where |Ψ 0 is the initial state of the many-body system at time t = 0, the post-quench HamiltonianĤ governs the time evolution for times t > 0, and we set = 1 from now on. The Loschmidt amplitude resembles the partition function of a thermal system with Hamiltonian H, however, the inverse temperature is replaced by the imaginary time τ = it, and an average is taken with respect to the initial state |Ψ 0 . In equilibrium settings, a thermal phase transition occurs, if a system is cooled below its critical temperature, and it abruptly changes from a disordered to an ordered phase. Similarly, dynamical phase transitions occur at critical times, when a quenched system suddenly changes from one phase to another with fundamentally different properties. Such transitions are manifested in the rate function which is the non-equilibrium analogue of the free energy density. In some cases, dynamical phase transitions occur, if a system is quenched across an underlying equilibrium phase transition, however, generally, there is no simple relation between dynamical and equilibrium phase transitions. In Fig. 1a, the system size, denoted by L, is the total number of spins along the chain. In the thermodynamic limit of infinitely large systems, dynamical phase transitions give rise to singularities in the rate function, for example a cusp as shown in Fig. 1b. However, this non-analytic behavior typically becomes apparent for very large systems, and it is hard to pinpoint for smaller systems. For this reason, dynamical phase transitions are difficult to capture in computations and simulations, where the numerical costs grow rapidly with system size.
Here we build on recent progress in Lee-Yang theories of thermal phase transitions [44][45][46][47] and use Loschmidt cumulants to predict dynamical phase transitions in strongly-correlated many-body systems using remarkably small system sizes. The Lee-Yang formalism of classical equilibrium phase transitions considers the zeros of the partition function in the complex plane of the external control parameters [48][49][50][51]. In a similar spirit, we treat the Loschmidt amplitude as a function of the complex-valued variable τ . The Loschmidt amplitude of a finite system is an entire function, which can be factorized as [52] where α is a constant, and τ k are the complex zeros of the Loschmidt amplitude. For a thermal system, the values of the inverse temperature for which the partition function vanishes are known as Fisher zeros [53]. We refer to the zeros of the Loschmidt amplitude as Loschmidt zeros. For a finite system, the zeros are isolated points in the complex plane. However, they grow denser as the system size is increased, and in the thermodynamic limit, they coalesce to form continuous lines and regions. Their intersections with the imaginary τ axis determine the real critical times, at which the rate function becomes non-analytic, and dynamical phase transitions occur [11]. As such, this phenomenology resembles the classical Lee-Yang theory of thermal phase transitions [48][49][50][51]. The central task is thus to determine the Loschmidt zeros. To this end, we define the Loschmidt moments and cumulants of the HamiltonianĤ of order n as and where τ is the basepoint, at which the moments and cumulants are evaluated. For τ = 0, the Loschmidt moments reduce to the moments of the Hamiltonian with respect to the initial state as Ĥ n 0 = Ψ 0 |Ĥ n |Ψ 0 . At finite times, the Loschmidt moments are Ĥ n τ = Ψ 0 |Ĥ n |Ψ(τ ) / Ψ 0 |Ψ(τ ) , where |Ψ(τ ) = e −τĤ |Ψ 0 is the time-evolved state. The cumulants can be obtained from the moments using the standard recursive formula For our purposes, it is now convenient to define the normalized Loschmidt cumulants κ n (τ ) as having used Eq. (3) to express them in terms of the zeros. This expression shows that the Loschmidt cumulants are dominated by the zeros that are closest to the (complex) basepoint τ , while the contributions from other zeros rapidly fall off with their inverse distance from the basepoint to the power of the cumulant order n. The main idea is now to extract the m closest zeros from 2m high Loschmidt cumulants, which we can calculate. For m = 2, this can be done by adapting the method from Refs. [44][45][46][47]. However, for arbitrary m, we use the general approach presented in App. A-B. For the systems that we consider in the following, we extract the m = 7 zeros closest to the movable basepoint using Loschmidt cumulants of order n = 9 to n = 22.
It should be emphasized that our approach hardly makes any assumptions about the quantum many-body system at hand or the method used for obtaining the cumulants. As outlined in App. C, we use a Krylov subspace method [54,55] to perform the complex time evolution and evaluate the Loschmidt moments and cumulants, which we then use for extracting the Loschmidt zeros. All results presented below have been obtained on a standard laptop, and the method can readily be adapted to a variety of interacting quantum many-body systems, also in higher dimensions.

III. INTERACTING KITAEV CHAIN
We first consider the spin-1/2 XYZ chain or, equivalently, the interacting Kitaev chain. The non-interacting limit maps to the XY model, which was solved exactly in the pioneering work of Ref. [3]. Here, we use Loschmidt cumulants to predict a dynamical quantum phase transition in the strongly interacting regime. The Hamiltonian of the spin-1/2 XYZ chain with a Zeeman field readŝ whereŜ α j are the spin-1/2 operators for the α = x, y, z component of the spin on site j of the chain of length L, the exchange couplings are denoted by J α , and h is the Zeeman field. We use twisted boundary conditions, andŜ z L+1 =Ŝ z 1 , where Φ is the twist angle. In the fermionic representation, obtained by a Jordan-Wigner transformation [56], the model maps to the interact-ing Kitaev chain of spinless fermions with operatorŝ c j andĉ † j , where the twist angle now enters in the last line through the parameter s Φ , which is +1, if the twist angle is Φ = 0, and −1, if Φ = π. These are the only two values of the twist angle used here. The parameters of the two Hamiltonians are related as J = −(J x + J y )/2, ∆ = (J y − J x )/2, µ = −h, and V = J z . Moreover, the number operator on site j isn j =ĉ † jĉ j , whileP = exp(iπ jn j ) is the parity operator.
The Kitaev chain describes a one-dimensional superconductor with a p-wave pairing term that is proportional to ∆, supporting two distinct topological phases.
The two values of the twist angle, Φ = 0, π, physically correspond to a magnetic flux equal to zero or half flux quantum threaded through the ring-shaped chain. These are the only distinct flux values that are consistent with superconducting flux quantization. It is useful to vary the boundary conditions, since in the non-interacting case (V = 0), which corresponds to the exactly solvable spin-1/2 XY model, the Loschmidt zeros can be labeled by the quasi-momentum k m = (2πm − Φ)/L with m = 0, . . . , L − 1 [3]. Thus, by using the two different values of Φ, we can sample the thermodynamic line of zeros twice as densely for a given system size. It turns out that even in the interacting case (V = 0) it is useful to vary the boundary conditions for the same reason.
We are now ready to investigate dynamical quantum phase transitions in the interacting Kitaev chain. To this end, we take for the initial state |Ψ 0 the ground state of the Hamiltonian (10) with |µ/J| > 1, which corresponds to the topologically trivial phase, and perform a quench into the topological regime with |µ/J| < 1 for later times, t > 0. As shown in Fig. 2, from the Loschmidt cumulants we can find the complex zeros of the Loschmidt amplitude even with attractive (V < 0) or repulsive (V > 0) interactions, for which an analytic solution is not available. In the left column, we first consider the non-interacting case, where the thermodynamic lines of zero can be determined analytically [3]. In panels a and b, we show the zeros found from the Loschmidt cumulants as the basepoint τ is moved along the paths denoted by A and B, respectively, while panel c shows the combined results. Remarkably, the Loschmidt cumulants allow us to map out the thermodynamic lines of zeros using chains of rather short lengths, L = 7 − 20, and thereby identify the crossing points with the imaginary axis, corresponding to the real critical times, where a dynamical phase transition occurs. The comparison between the exact and the approximate zeros obtained from the Loschmidt cumulants provides an important estimate of the accuracy. In the worst cases, the accuracy is an order of magnitude better than the size of the markers in Fig. 2 (see App. B). We note that our choice of the paths in Fig. 2 was guided by our knowledge of the zeros in the non-interacting case. However, more generally, without any specific knowledge of a system, one may choose paths that scan the complex plane, in particular along the imaginary time axis and its immediate vicinity (since those zeros determine whether and when the system exhibits a dynamical phase transition).
Having benchmarked our approach in the noninteracting case, we move on to the strongly interacting regime. In the second column of Fig. 2, we show the Loschmidt zeros for repulsive interactions, which tend to shift the critical crossing point with the imaginary axis to earlier times. A more dramatic effect is observed in the third and fourth columns, where we gradually increase the attractive interactions. In this case, the dynamical phase transition happens at later times, and eventually, for sufficiently strong interactions, the thermodynamic line of zero no longer crosses the imaginary axis, implying the absence of a dynamical phase transition.
While in the noninteracting limit the small systems reproduce the thermodynamic lines essentially exactly, interactions give rise to finite-size effects when two different lines come close. Despite that, sufficiently isolated lines and segments, such as the ones that determine the dynamical phase transitions in Fig. 2, remain scale-invariant. We stress that these results are obtained for very small chains of lengths from L = 10 to L = 20, which, while remarkable, is in line with similar observations for the Lee-Yang zeros in classical equilibrium systems [44][45][46][47]. In particular, for strongly interacting systems, the use of such system sizes makes the approach very attractive from a computational point of view, since direct calculations of the Loschmidt amplitude typically require system sizes that are an order of magnitude larger, in generic cases with an exponential increase in the computational cost.

IV. SPIN-1 HEISENBERG CHAIN
The Kitaev chain from above possesses an exactly solvable limit, which provides an important benchmark for the use of the Loschmidt cumulants. However, generically, exact solutions are not available, which makes the usefulness of the Loschmidt cumulants further evident. For this reason, we now consider the spin-1 Heisenberg chain, which harbors a rich phase diagrams both with symmetry-broken phases and a topological phase, the Haldane phase [58]. The spin-1 Heisenberg chain is defined by the Hamiltonian whereŜ i are spin-1 operators, the exchange couplings between neighboring spins are denoted by J and J z , while D and B characterize the single-spin uniaxial anisotropy and the biquadratic exchange coupling, respectively. The first line defines the spin-1 XXZ model, while for J = J z = 3B and D = 0, one obtains the Affleck-Kennedy-Lieb-Tasaki (AKLT) model, whose ground state is known In panels b1, b2, c1, c2 the crossings of the effective thermodynamic lines of zeros with the imaginary axis are shown and compared again with the critical times obtained in Ref. [27].
explicitly [57], despite the fact that the Hamiltonian (11) is not exactly solvable. Again, we employ twisted boundary conditions as defined in Eq. (9) for five different values of the phase Φ = 0, π/4, π/2, 3π/4, π. We consider two kinds of quenches in which different parameters in the Hamiltonian (11) are abruptly changed at t = 0.
In the first quench, we initialize the system in the AKLT ground state, which is a representative of the topological Haldane phase, and evolve it with the Hamiltonian (11)  . a, Loschmidt zeros for the quench to Jz = 1, which is not sufficiently large to reach the Néel phase. In this case, the zeros do not cross the imaginary axis, and no dynamical phase transition occurs. b, c, d Similar to panel a, but with Jz = 2, 3, 4. In this case, several dynamical phase transitions occur as shown for example in panels e1-3 and f -g. The critical times are shown as red crosses and are estimated as done in Fig. 3 using the zeros for Φ = π/2 only (see App. D for details).
phase. The same quenches have been explored in Ref. [27] for system sizes up to L = 120 using matrix product states, providing us with an important benchmark. Figure 3 shows the Loschmidt zeros for finite system sizes extracted from the Loschmidt cumulants for the first quench. We use twisted boundary conditions to gauge finite-size effects as the position of the Loschmidt zeros is expected to become insensitive to the phase Φ for very large systems. By contrast, for the relatively small system sizes used in Fig. 3, finite-size effects are pronounced in particular in panel a, which shows the zeros for the D/J = 2 quench. This value is the closest to the critical one D c /J 1 (with J = J z and B = 0) separating the Haldane phase from the large-D phase [58], provid-ing a plausible reason for the enhanced finite-size effects. Importantly, as discussed in App. D, the oscillatory pattern of zeros for different system sizes and twist angles is highly regular, which enables us to filter out the finitesize effects. In this prescription, a thermodynamic line of zeros is approximated by the smooth line of zeros emerging at the twist angle Φ = π/2.
The critical times of the transition, obtained from the crossings of the thermodynamic lines of zeros with the imaginary axis (see panels a1, b1-2 and c1-2 in Fig. 3), are in excellent agreement with the critical times obtained directly from the Loschmidt amplitude that was calculated using state-of-the-art computations in Ref. [27]. However, in contrast to Ref. [27], which con-siders nearly an order of magnitude larger systems, our results are obtained from chain lengths up to L = 16. This comparison provides an illustration of the power of our method in treating strongly correlated many-body systems.
While the exact correspondence between dynamical phase transitions and the equilibrium phase transitions of the respective model remains unknown [11], dynamical phase transitions are often observed, when the ground states of the initial and final Hamiltonians belong to different equilibrium phases. To explore this general scenario in the case of transitions between a topological phase and a symmetry-broken phase in a strongly correlated system, we solve for the first time quenches between the topological Haldane phase and the symmetry-broken Néel phase [58]. In Fig. 4, we depict the Loschmidt zeros for the initial state with D = B = 0 and quenching J z from J z /J = 1/2 to the final values J z /J = 1, 2, 3, 4. The equilibrium quantum phase transition occurs at the critical value J z,c 1.2J [58]. Indeed, our results confirm that no dynamical phase transition is observed when J z /J = 1, since all the Loschmidt zeros have negative real part as shown in panel a of Fig. 4. By contrast, for the other final values of J z , which would put the equilibrium system in the antiferromagnetic Néel phase, dynamical phase transitions are observed. As in the first quench, finite-size effects are suppressed for quenches, where the final state resides deeper in the gapped phase.
As we see in panels e1-3 of Fig. 4, the Loschmidt zeros for different system sizes and boundary conditions have a structure similar to the one observed in the Haldaneto-large-D quench. The same prescription as above irons out the finite-size oscillations and results in a smooth approximation of the thermodynamic lines of zeros. The critical times can then be accurately read off from the data obtained for chain lengths of L ≤ 16, as in the case of the first quench.

V. EXPERIMENTAL PERSPECTIVES
In the previous sections, we focused on using the Loschmidt cumulants for predicting dynamical phase transitions based on numerical calculations. However, as we will now discuss, our approach also provides perspectives for future experiments. We will show that it is possible to predict the first critical time of a quantum many-body system by measuring the fluctuations of the energy in the initial state. We will also discuss the prospects of implementing our method on a near-term quantum computer with a small number of qubits.
The Loschmidt moments are generally complex-valued, and it is not obvious how they can be measured. However, at the initial time, τ = 0, the Loschmidt moments simplify to the moments of the post-quench Hamiltonian with respect to the initial state as Ĥ n 0 = Ψ 0 |Ĥ n |Ψ 0 . Thus, by repeatedly preparing the system in the state |Ψ 0 and measuring the energy given by the post-quench HamiltonianĤ, one can construct the distribution of the energy and extract the corresponding moments and cumulants. From the cumulants, it is then possible to extract the closest Loschmidt zeros as demonstrated in Fig. 5 following a quench in a Heisenberg chain of lengths L = 5, . . . , 9. From these results, we predict the critical time to be around t c 0.42 as indicated by a red cross. This perspective is fascinating: by measuring the initial energy fluctuations, it is possible to predict the later time at which a dynamical phase transition will occur. The idea behind such an experiment does not depend in detail on the actual physical implementation, and from a practical point of view different platforms may provide certain advantages. We expect, for example, that an experiment could be realized with atoms in optical lattices [59] or with spin chains on surfaces [60], systems that both offer a high degree of control and flexibility. As illustrated in Fig. 5, it will be necessary to measure the high cumulants of the energy fluctuations. For large systems, accurate measurements of high cumulants are challenging, since the central-limit theorem dictates that distributions tend to be Gaussian with nearly vanishing high cumulants. However, for the small systems that we consider, the situation is different, and several quantum transport experiments have measured cumulants of up to order 20 [61,62] and used them for determining the zeros of generating functions [63,64], which are similar to the Loschmidt amplitude. Thus, an experimental determination of Loschmidt zeros for small interacting quantum systems appears feasible with current technology.
Our method may also be implemented on small nearterm quantum computers, which are now becoming available. Such quantum computers allow for the specific tailoring of any desired Hamiltonian and for timeevolving an initial state both in real and imaginary time [65,66]. Thus, it will be possible to evaluate a time-evolved state of the form |Ψ(τ ) = e −τĤ |Ψ 0 and subsequently calculate the Loschmidt moments Ĥ n τ = Ψ 0 |Ĥ n |Ψ(τ ) / Ψ 0 |Ψ(τ ) and the corresponding cumulants from which the Loschmidt zeros are obtained. Again, the favorable scaling properties of our method become important, as they make it possible to predict the critical times of a quantum many-body system with only 10 to 20 constituents. Such sizes can soon be simulated on quantum computers with a limited number of qubits.

VI. CONCLUSIONS
We have demonstrated that Loschmidt cumulants are a powerful tool to unravel dynamical phase transitions in strongly interacting quantum many-body systems after a quench, making it possible to accurately predict the critical times of a quantum many-body system using remarkably small system sizes. Using modest computational power, we have explored dynamical phase transitions in the Kitaev chain and the spin-1 Heisenberg chain with a specific focus on the role of strong interactions. As we have shown, our approach circumvents the existing bottleneck of computing the full non-equilibrium dynamics of large quantum many-body systems, and instead we track the zeros of the Loschmidt amplitude in the complex plane of time in a similar spirit to the classical Lee-Yang theory of equilibrium phase transitions. As such, our approach paves the way for systematic investigations of the far-from-equilibrium properties of interacting quantum many-body systems, and we foresee many exciting perspectives ahead. In particular, our method can immediately be applied to dynamical phase transitions in dimensions higher than one, and the ease of implementing it may be critical for comprehensive investigations of the finite-size scaling close to a dynamical phase transition. We have also shown that our approach paves the way for exciting experimental developments by making it possible to predict the first critical time of a quantum many-body system in the thermodynamic limit by measuring the initial energy fluctuations in a much smaller system. In addition, due to the favorable scaling of our method, it seems feasible that it can be implemented on a near-term quantum computer with a limited number of qubits. In a broader perspective, the advances presented here may not only be useful for understanding the dynamical non-equilibrium properties of large quantum systems. They may also be helpful in designing novel quantum materials with specific, desired properties.

ACKNOWLEDGMENTS
We thank the authors of Ref. [27] for providing us with their results for the spin-1 Heisenberg chain, which we used to extract the critical times indicated in Fig. 3. Here we present the basic idea of the method used to extract the Loschmidt zeros from the Loschmidt cumulants. More details on the specific procedure employed to obtain the results shown in Figs. 2-4 are provided in App. B, where we for instance explain how we estimate the error that affects the approximate zeros extracted with our method.
In Eqs. (3,7), repeated zeros are allowed such that every distinct zero appears in the series (7) as many times as its multiplicity. It is convenient in the following to use a different labelling of the Loschmidt amplitude zeros in which the index k runs over the distinct zeros (τ k = τ k for k = k ), and also to denote by d k the multiplicity of the k-th zero. Using this convention, Eq. (7) becomes On the right hand side, we have introduced λ k = 1/(τ k − τ ) and truncated the sum to the m zeros closest to the basepoint. This is a good approximation for large n since the contribution of each zero to the normalized cumulant κ n (τ ) is suppressed by its inverse distance to the basepoint raised to the power of the cumulant order. We now determine the zeros of the Loschmidt amplitude based on the fact that, if Eq. (A1) were exact, the normalized cumulants would satisfy a homogeneous linear difference equation of degree m of the form κ n = a 1 κ n−1 + a 2 κ n−2 + . . . + a m (A2) for some coefficients a l . Indeed, the general solution of Eq. (A2) is given by the right hand side of Eq. (A1), where λ k are the characteristic roots of the characteristic equation associated to Eq. (A2), see Eq. (A4) below, and d k can be arbitrary coefficients due to linearity. This observation is crucial for inverting Eq. (A1) and extracting the zeros from the cumulants as we explain below. We note that Eqs. (A1,A2) are exact only if the Loschmidt amplitude Z(τ ) is a polynomial with m distinct zeros. In this case, the method provides exactly all the zeros of Z(τ ), independently of the cumulant orders used. In the general case, where Z(τ ) is an entire function, the method provides approximate zeros τ (n) k that depend on the cumulant orders and converge to the exact zeros for n → ∞. Associated to the approximate zeros, the method provides also a sequence of approximate multiplicities d (n) k for each k, which converges to the exact multiplicity d k .
The first step of the method is to compute the coefficients a l=1,...,m in Eq. (A2). This can be done by solving a linear system of m equations, which requires the knowledge of 2m consecutive cumulants (κ l , with n − m ≤ l ≤ n + m − 1) and takes the form The square matrix on the left hand side is a Toeplitz matrix. With the notation a which is a polynomial equation in λ. This provides the m characteristic roots λ where τ are two sequences that converge for n → ∞ to the two closest zeros τ 0 and τ 1 ( = τ * 0 in general). The third and final step is the calculation of the coefficients d where we have dropped the superscript from λ k → d k for n → ∞ together with the approximate zeros τ (n) k → τ k . We use this fact to select the approximate zeros which are the best approximations of the exact zeros when the basepoint is varied as discussed in App. B.

Appendix B: Extracting Loschmidt zerosnumerical convergence and error estimates
In order to obtain accurate approximations of the exact zeros, one can increase the cumulant order (the parameter n in Eq. (A3)) to observe the convergence of the approximate zeros τ (n) k . However, we use a different procedure in our work, which is almost automatic and has proven particularly effective for investigations of dynamical quantum phase transitions. The idea is based on the fact that the approximate multiplicities obtained from Eq. (A6) converge to the exact multiplicities d k − | < r, where r is a fixed threshold and is a chosen integer. This condition appears to be necessary but not sufficient to guarantee that τ (n) k is a good approximation of τ k . Indeed, it can occur that the approximate zero of a pair (τ is not a good approximation of any exact zero even if the condition |d (n) k − | < r is satisfied for some integer . However, this kind of false positives are in practice quite rare for small enough r and can be easily detected and discarded by applying the method at different basepoints, as explained below. The variant of the method in the case of a pair of conjugate zeros (τ 0 = τ * 1 , m = 2) has been applied earlier to study critical phenomena in classical equilibrium problems [44][45][46][47], and here we have extended the method to an arbitrary number of zeros, which are not necessarily pairwise conjugate. These are crucial advancements that have allowed us to apply the cumulant method to the study of dynamical quantum phase transitions. The use of the coefficients d (n) k for selecting the best approximate zeros is also a new technique, which makes the method very practical and efficient for the prediction of dynamical quantum phase transitions and is used for the first time in this work.
The results shown in Figs. 2-4 have been obtained by evolving the initial state along the imaginary axis in the complex τ -plane (path A in Fig. 2a). In the case of the Kitaev chain, we perform the time evolution also along path B, see Fig. 2b, a straight path starting at τ = 0 and ending at τ = −4 + 6i. In this way we can map out zeros in different regions of the complex plane. The intermediate values of τ along the evolution path are the basepoints at which the cumulants are computed. We use a fine grid in which adjacent basepoints are at a distance δτ = 0.001. For each basepoint τ we apply Eqs. (A3), (A4) and (A6) using cumulants κ n (τ ) from n = 9 to n = 22, and we obtain m = 7 distinct approximate zeros τ (j) k=0,...,6 and corresponding coefficients d k=0.,...,6 . Here we change the notation slightly: the superscript in the pair (τ k ) labels the distinct basepoints here. This is because the cumulant orders used in the method are kept fixed in this case and only the basepoint is varied along the path.
As explained above, we select only the pairs (τ k − 1| < r = 0.01. We have not found zeros with multiplicity > 1 in the systems considered in our work. After this selection step, one can visually verify that the approximate zeros tend to agglomerate in well separated clusters. In the case where the Loschmidt amplitude can be computed analytically and the exact zeros are known (the Kitaev chain with V = 0 [3], Fig. 2a-c), one can also see that each cluster corresponds to one of the exact zeros. We use the standard k-means++ clustering algorithm [67] to classify the approximate zeros in distinct clusters. This algorithm requires to introduce manually the number of clusters, which can be estimated by visual inspection. Each cluster obtained in this way typically consists of hundreds or thousands of pairs (τ (j) k , d (j) k ). Clusters with less then ten pairs are discarded to eliminate any false positives.
Finally, within each cluster we select the pair (τ k −1| takes its smallest value. The approximate zero τ (j) k of this pair gives the best estimate of the location of one zero. Indeed, we have observed in exactly solvable cases that the same pair minimizes also the distance |τ (j) k − τ k | from the closest exact zero. Notice that one can generally resolve more than m zeros in a single run of time evolution by using the above procedure since different zeros are resolved for different basepoints.
The standard deviation s of the real and imaginary parts of the approximate zeros within each cluster provides a rough estimate of the error, i.e. the distance from the exact zero. Typically, we obtain s ≈ 10 −3 using the procedure presented above. By comparison with the exact solution we have verified that this is a reasonable estimate or even an overestimation in most cases. Another way to estimate the error is to compare the zeros obtained by evolving along different paths as in Fig. 2ac. Often the same exact zero can be resolved by using both paths and thus one obtains two different estimates of its location whose distance is an estimate of the er-ror. The error estimate obtained in this way turns out to be essentially the same as the one obtained from the cluster standard deviation. An error of order 10 −3 is not visible on the scale of Figs. 2-4, since it is an order of magnitude smaller than the size of the markers. Therefore, the fact that in the interacting case the zeros seem not to organize in well defined lines (in contrast to the non-interacting case in Fig. 2a-c) has to be entirely attributed to finite-size effects and not to the approximate nature of the zeros obtained with our method.

Appendix C: Krylov subspace method
The evolution along a straight path in the complex τ -plane is performed by using the standard Krylov subspace method [54,55]. A time step δτ in the evolution is performed by first computing an orthonormal basis B = {|v 0 = |Ψ(τ ) , |v 1 , . . . , |v Nvec } of the Krylov subspace We use the QR decomposition to compute the orthonormal basis of K Nvec to ensure that there is no loss of orthogonality as in the standard Lanczos algorithm [55]. Then, the approximate time-evolved state is obtained as whereĤ eff is the effective Hamiltonian, that is an operator which acts on the Krylov subspace and whose matrix elements are given by v i |Ĥ |v j with |v i ∈ B. It is represented by a square matrix of dimension N vec + 1 whose exponential is easy to evaluate since we take N vec = 8 in our case. As suggested in Section 5 of Ref. [55], the effective Hamiltonian is forcibly set to be a tridiagonal matrix to improve stability. In order to estimate the accuracy of Eq. (C2) we perform the time evolution on two Kyrlov subspaces K N vec and K N vec with N vec + 1 = N vec ≤ N vec and compute the distance between the approximate evolved states d = |Ψ (τ + δτ ) − |Ψ (τ + δτ ) 2 . If d < 10 −10 the state |Ψ (τ + δτ ) is stored and used to evaluate the moments of the Loschmidt amplitude according to Ĥ n τ +δτ = Ψ 0 |Ĥ n |Ψ (τ + δτ ) . On the other hand, if for a given δτ the condition is not satisfied even for N vec = N vec , the Krylov subspace K Nvec in Eq. (C1) is recomputed using as the seed state the last stored vector |Ψ (τ + δτ ) . The whole process is repeated up to the desired final τ .
An important advantage of the Krylov time evolution algorithm described above is that intermediate values of τ along the evolution path are easily accessible at a negligible computational cost. We take advantage of this fact to compute the cumulants on a fine grid of basepoints along the evolution path. The computationally expensive part of the method is the calculation of the Krylov subspace (C1). In order to compute the moments and the cumulants, a larger Krylov subspace (N vec = 22 in our case) has to be calculated only at τ = 0. This is a small numerical overhead since the Krylov subspace has to be recalculated many times during the evolution. With the specific parameters given above, we find by comparison with exactly solvable cases that the cumulants are computed with a relative precision of 10 −9 , which is sufficient for our purposes.

Appendix D: Twisted boundary conditions
The twisted boundary conditions that we employ are known as a tool to filter out finite-size effects. For example, twisted boundary conditions have been employed to analyze energy-level crossings related to quantum phase transitions in spin-1 systems [68]. Here we discuss how these boundary conditions can also help to gauge finitesize effects in the Loschmidt zeros. For the Haldane chain, Eq. (11), the spectrum exhibits a 1/L dependence on the boundary conditions [68]. Thus, the Loschmidt zeros, determined by the spectrum of the Hamiltonian, are independent of the boundary conditions for sufficiently large system sizes. This is a reasonable assumption for most naturally occurring systems. However, one may in principle envision situations that deviate from this generic behavior, and one may even construct specific examples that do [69]. On the other hand, in finite systems, the many-body energies E i (Φ) depend on the twist angle Φ, in particular the highly excited states. The Loschmidt zeros inherit this property -a single zero, say τ k , gives rise to a multiplet of zeros τ k (Φ) corresponding to different values of Φ as seen in Fig. 3 (in particular, see panel a1) and in Fig. 4. Fortunately, we observe a regular pattern, which can be employed as a diagnostic tool. Specifically, all E i (Φ) are even functions of Φ, having their extreme values at Φ = 0 or Φ = π. As seen in Fig. 3a1, the boundary points of the multiplets of zeros exactly correspond to Φ = 0 or Φ = π. In the thermodynamic limit, these multiplets converge to a single point, τ k , for all angles Φ.
A natural question is how one can obtain the best approximation for the thermodynamic line of zeros. It is reasonable to expect that the line is situated between the extreme zeros at the twist angles Φ = 0 and π. The curve of zeros corresponding to a general value of Φ exhibits size-dependent oscillations around the mean value as seen in Fig. 3a1 by the zig-zag line corresponding to Φ = 0. However, these finite-size oscillations are suppressed for zeros corresponding to the twist angle Φ = π/2, which form a smooth curve. This observation suggests that the twist angle Φ = π/2 is special, and that it is the best approximation for the thermodynamic line of zeros. This prescription is independently supported by the remarkable agreement between the critical times that we obtain and the results of Ref. [27], which are also indicated in Fig. 3. A similar pattern is observed for the quench from