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

Sebastiano Peotta  Computational Physics Laboratory, Physics Unit, Faculty of
Engineering and Natural Sciences, Tampere University, FI-33014 Tampere,
Finland  Department of Applied Physics, Aalto University, FI-00076 Aalto,
Finland  Helsinki Institute of Physics, FI-00014 University of Helsinki,
Finland  Fredrik Brange  Department of Applied Physics, Aalto University,
FI-00076 Aalto, Finland  Aydin Deger  Department of Applied Physics, Aalto
University, FI-00076 Aalto, Finland  Teemu Ojanen  Email<EMAIL_ADDRESS>Computational Physics Laboratory, Physics Unit, Faculty of Engineering and
Natural Sciences, Tampere University, FI-33014 Tampere, Finland  Helsinki
Institute of Physics, FI-00014 University of Helsinki, Finland  Christian
Flindt  Department of Applied Physics, Aalto University, FI-00076 Aalto,
Finland

###### Abstract

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 Chandler (1987); Goldenfeld (2018), 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 Heyl
_et al._ (2013); Karrasch and Schuricht (2013); Andraschko and Sirker (2014);
Vajna and Dóra (2014, 2015); Budich and Heyl (2016); Zvyagin (2016); Halimeh
and Zauner-Stauber (2017); Heyl (2018). 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 one-dimensional 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 trapped ions Jurcevic  _et al._ (2017); Zhang
_et al._ (2017), optical lattices Fläschner  _et al._ (2018), and several
other systems that offer a high degree of control Guo  _et al._ (2019); Tian
_et al._ (2019); Wang  _et al._ (2019); Tian  _et al._ (2020); Xu  _et al._
(2020). 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\simeq 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 Pozsgay (2013);
Karrasch and Schuricht (2013); Andraschko and Sirker (2014); Kriel  _et al._
(2014); Sharma  _et al._ (2015); Halimeh and Zauner-Stauber (2017); Zauner-
Stauber and Halimeh (2017); Homrighausen  _et al._ (2017); Heyl  _et al._
(2018); Kennes  _et al._ (2018); Hagymási  _et al._ (2019); Lacki and Heyl
(2019); Huang  _et al._ (2019); Halimeh  _et al._ (2020); Žunkovič  _et al._
(2018); Weidinger  _et al._ (2017); Feldmeier  _et al._ (2019). 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 Heyl  _et al._ (2013); Karrasch and Schuricht (2013); Andraschko and
Sirker (2014); Vajna and Dóra (2014); Heyl (2015); Vajna and Dóra (2015);
Schmitt and Kehrein (2015); Budich and Heyl (2016); Zvyagin (2016); Halimeh
and Zauner-Stauber (2017); Fogarty  _et al._ (2017); Bhattacharya and Dutta
(2017); Heyl (2018); Kennes  _et al._ (2018); Sedlmayr  _et al._ (2018);
Jafari  _et al._ (2019); Defenu  _et al._ (2019); Najafi  _et al._ (2019);
Gulácsi  _et al._ (2020); Zamani  _et al._ (2020). 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 related 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.

Figure 1:  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,
$t_{c,1},t_{c,2}\ldots$, 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 $\tau$.

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\simeq 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,

$\mathcal{Z}(it)=\left\langle\Psi_{0}\right|e^{-it\mathcal{\hat{H}}}\left|\Psi_{0}\right\rangle,$
(1)

where $\left|\Psi_{0}\right\rangle$ is the initial state of the many-body
system at time $t=0$, the post-quench Hamiltonian $\mathcal{\hat{H}}$ governs
the time evolution for times $t>0$, and we set $\hbar=1$ from now on. The
Loschmidt amplitude resembles the partition function of a thermal system with
Hamiltonian $\mathcal{\hat{H}}$, however, the inverse temperature is replaced
by the imaginary time $\tau=it$, and an average is taken with respect to the
initial state $\left|\Psi_{0}\right\rangle$. 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

$\lambda(t)=-\frac{1}{L}\ln|\mathcal{Z}(it)|^{2},$ (2)

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 Deger  _et al._ (2018); Deger and Flindt (2019, 2020); Deger  _et
al._ (2020) 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 Yang and Lee (1952); Lee and Yang (1952); Blythe and Evans
(2003); Bena  _et al._ (2005). In a similar spirit, we treat the Loschmidt
amplitude as a function of the complex-valued variable $\tau$. The Loschmidt
amplitude of a finite system is an entire function, which can be factorized as
Arfken  _et al._ (2012)

$\mathcal{Z}(\tau)=e^{\alpha\tau}\prod_{k}\left(1-\tau/\tau_{k}\right),$ (3)

where $\alpha$ is a constant, and $\tau_{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 Fisher . 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 $\tau$ axis determine the real critical
times, at which the rate function becomes non-analytic, and dynamical phase
transitions occur Heyl (2018). As such, this phenomenology resembles the
classical Lee-Yang theory of thermal phase transitions Yang and Lee (1952);
Lee and Yang (1952); Blythe and Evans (2003); Bena  _et al._ (2005).

The central task is thus to determine the Loschmidt zeros. To this end, we
define the Loschmidt moments and cumulants of the Hamiltonian
$\mathcal{\hat{H}}$ of order $n$ as

$\langle\mathcal{\hat{H}}^{n}\rangle_{\tau}=(-1)^{n}\frac{\partial_{\tau}^{n}\mathcal{Z}(\tau)}{\mathcal{Z}(\tau)}$
(4)

and

$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{\hat{H}}^{n}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\tau}=(-1)^{n}\partial_{\tau}^{n}\ln\mathcal{Z}(\tau),$
(5)

where $\tau$ is the basepoint, at which the moments and cumulants are
evaluated. For $\tau=0$, the Loschmidt moments reduce to the moments of the
Hamiltonian with respect to the initial state as
$\langle\mathcal{\hat{H}}^{n}\rangle_{0}=\left\langle\Psi_{0}\right|\mathcal{\hat{H}}^{n}\left|\Psi_{0}\right\rangle$.
At finite times, the Loschmidt moments are
$\langle\mathcal{\hat{H}}^{n}\rangle_{\tau}=\left\langle\Psi_{0}\right|\mathcal{\hat{H}}^{n}\left|\Psi(\tau)\right\rangle/\langle\Psi_{0}|\Psi(\tau)\rangle$,
where
$\left|\Psi(\tau)\right\rangle=e^{-\tau\mathcal{\hat{H}}}\left|\Psi_{0}\right\rangle$
is the time-evolved state. The cumulants can be obtained from the moments
using the standard recursive formula

$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{\hat{H}}^{n}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\tau}=\langle\mathcal{\hat{H}}^{n}\rangle_{\tau}-\sum_{m=1}^{n-1}\binom{n-1}{m-1}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{\hat{H}}^{m}\mathclose{\hbox{\set<EMAIL_ADDRESS>(6)

For our purposes, it is now convenient to define the normalized Loschmidt
cumulants $\kappa_{n}(\tau)$ as

$\kappa_{n}(\tau)=\frac{(-1)^{n-1}}{(n-1)!}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{\hat{H}}^{n}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\tau}=\sum_{k}\frac{1}{(\tau_{k}-\tau)^{n}},\quad
n>1,$ (7)

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 $\tau$, 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. Deger  _et al._ (2018);
Deger and Flindt (2019, 2020); Deger  _et al._ (2020). 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 Park and
Light (1986); Paeckel  _et al._ (2019) 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.

Figure 2: The interacting Kitaev chain. We quench the chemical potential from
the trivial phase $\mu=-1.4$ to the topological phase $\mu=-0.2$ for $t>0$
with fixed $\Delta=0.3$ (all parameters and the inverse time $\tau^{-1}$ are
expressed in units of $J=1$). a-c, Complex zeros for different system sizes
($L=7-20$) and boundary conditions ($\Phi=0,\pi$) in the noninteracting case,
compared with the exact thermodynamic lines of zeros. Only the zeros within a
finite range from the basepoint can be accurately obtained from the Loschmidt
cumulants. This fact is illustrated by moving the basepoint $\tau$ along two
different paths (paths A and B in panels a and b). Panel c combines the
results from panels a and b. d-f, Loschmidt zeros and critical times obtained
with repulsive interactions ($V>0$) along the two paths. The lines of zeros
for $V=0$ (dashed line) are shown for comparison. The critical time $t_{c}$,
shown in each panel as a red cross, is obtained as the intersection between
the imaginary axis and the line drawn from the zero $\tau_{-}$ with smallest
negative real part (in absolute value) to the zero $\tau_{+}$ with the
smallest positive real part. The error on the critical time is estimated as
$\Delta
t_{c}=\max(|t_{c}-\mathrm{Im}\,\tau_{-}|,|t_{c}-\mathrm{Im}\,\tau_{+}|)$. g-l,
Evolution of zeros and critical times with increasing attractive interactions
($V<0$). For very strong interactions ($V=-1$, panel l), the zeros do not
cross the imaginary axis, signalling the absence of a dynamical quantum phase
transition. As discussed in App. B, the numerical error in the zeros is of the
order of $10^{-3}$.

##  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. Heyl  _et al._ (2013). 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 reads

$\mathcal{\hat{H}}=\sum_{\alpha,j=1}^{L}J_{\alpha}\hat{S}_{j}^{\alpha}\hat{S}_{j+1}^{\alpha}-h\sum_{j=1}^{L}\hat{S_{j}^{z}}\,,$
(8)

where $\hat{S}_{j}^{\alpha}$ are the spin-1/2 operators for the $\alpha=x,y,z$
component of the spin on site $j$ of the chain of length $L$, the exchange
couplings are denoted by $J_{\alpha}$, and $h$ is the Zeeman field. We use
twisted boundary conditions,

$\begin{split}\hat{S}_{L+1}^{x}&=\cos(\Phi)\hat{S}_{1}^{x}+\sin(\Phi)\hat{S}^{y}_{1},\\\
\hat{S}_{L+1}^{y}&=-\sin(\Phi)\hat{S}_{1}^{x}+\cos(\Phi)\hat{S}^{y}_{1},\end{split}$
(9)

and $\hat{S}_{L+1}^{z}=\hat{S}_{1}^{z}$, where $\Phi$ is the twist angle. In
the fermionic representation, obtained by a Jordan-Wigner transformation
Franchini (2016), the model maps to the interacting Kitaev chain of spinless
fermions with operators $\hat{c}_{j}$ and $\hat{c}_{j}^{\dagger}$,

$\mathcal{\hat{H}}=-\frac{1}{2}\sum_{j=1}^{L-1}\left(J\,\hat{c}_{j}^{\dagger}\hat{c}_{j+1}+\Delta\,\hat{c}_{j}^{\dagger}\hat{c}_{j+1}^{\dagger}+\text{H.c.}\right)+V\sum_{j=1}^{L-1}\left(\hat{n}_{j}-\frac{1}{2}\right)\left(\hat{n}_{j+1}-\frac{1}{2}\right)-\mu\sum_{j=1}^{L}\hat{c}_{j}^{\dagger}\hat{c}_{j}+\frac{s_{\Phi}\hat{P}}{2}\left(J\,\hat{c}_{L}^{\dagger}\hat{c}_{1}+\Delta\,\hat{c}_{L}^{\dagger}\hat{c}_{1}^{\dagger}+\text{H.c.}\right),$
(10)

where the twist angle now enters in the last line through the parameter
$s_{\Phi}$, which is $+1$, if the twist angle is $\Phi=0$, and $-1$, if
$\Phi=\pi$. 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$,
$\Delta=(J_{y}-J_{x})/2$, $\mu=-h$, and $V=J_{z}$. Moreover, the number
operator on site $j$ is $\hat{n}_{j}=\hat{c}_{j}^{\dagger}\hat{c}_{j}$, while
$\hat{P}=\exp(i\pi\sum_{j}\hat{n}_{j})$ is the parity operator.

The Kitaev chain describes a one-dimensional superconductor with a $p$-wave
pairing term that is proportional to $\Delta$, supporting two distinct
topological phases. The two values of the twist angle, $\Phi=0,\pi$,
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\pi m-\Phi)/L$
with $m=0,\dots,L-1$ Heyl  _et al._ (2013). Thus, by using the two different
values of $\Phi$, 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\neq 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
$\left|\Psi_{0}\right\rangle$ the ground state of the Hamiltonian (10) with
$|\mu/J|>1$, which corresponds to the topologically trivial phase, and perform
a quench into the topological regime with $|\mu/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 Heyl  _et al._ (2013). In panels
a and b, we show the zeros found from the Loschmidt cumulants as the basepoint
$\tau$ 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 non-interacting 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 Deger  _et al._ (2018); Deger and
Flindt (2019, 2020); Deger  _et al._ (2020). 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.

Figure 3: The Heisenberg chain. We quench the system from the Haldane phase to
the large-$D$ phase. The initial state is the ground state of the model (11)
with $J=J_{z}=3B$ and $D=0$ (the AKLT state Affleck  _et al._ (1987)). For the
post-quench Hamiltonian, we set $B=0$, while $D$ can take the values 2, 3, or
4. Here $J=J_{z}=1$ is the unit of energy and inverse time. a, Loschmidt zeros
for $D=2$. Panel a1 is a magnified view of the area within the black rectangle
in panel a. From panels a and a1 one can clearly see how the position of a
Loschmidt zero for fixed $L$ depends on the twist angle $\Phi$, which is a
finite-size effect. It is also useful to consider a fixed twist angle and vary
the system size as in the case of the zeros connected by the dash-dotted line
in panel a1 ($\Phi=0$, $L=13$, $14$, $15$, $16$). The finite-size dependence
is suppressed for the zeros corresponding to the twist angle $\Phi=\pi/2$,
defining the effective thermodynamic line of zeros (solid line, see App. D).
The critical time, determined by the crossing of the effective line with the
imaginary axis, is in excellent agreement with the result of Ref. Hagymási
_et al._ (2019) (red bar) obtained using matrix product states (MPS). b, c
Same as in panel a but with $D=3,4$. Finite-size effects are suppressed with
increasing $D$. 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. Hagymási  _et al._ (2019).

##  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 Chen
_et al._ (2003). The spin-1 Heisenberg chain is defined by the Hamiltonian

$\begin{split}\mathcal{\hat{H}}=&\sum_{j=1}^{L}\left[J\big{(}\hat{S}_{j}^{x}\hat{S}_{j+1}^{x}+\hat{S}_{j}^{y}\hat{S}_{j+1}^{y}\big{)}+J_{z}\hat{S}_{j}^{z}\hat{S}_{j+1}^{z}\right]\\\
&+D\sum_{j=1}^{L}(\hat{S_{j}^{z}})^{2}+B\sum_{j=1}^{L}(\hat{\mathbf{S}}_{j}\cdot\hat{\mathbf{S}}_{j+1})^{2}\,,\end{split}$
(11)

where $\hat{S}^{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 explicitly Affleck  _et al._ (1987),
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 $\Phi=0,\,\pi/4,\,\pi/2,\,3\pi/4,\,\pi$. We consider two
kinds of quenches in which different parameters in the Hamiltonian (11) are
abruptly changed at $t=0$.

Figure 4: Quench from the Haldane phase to the Néel phase. The initial state
is the ground state of the model (11) with $J_{z}/J=1/2$ and $D=B=0$. The
quench is performed by changing the parameter $J_{z}$ to the values
$J_{z}=1,\,2,\,3,\,4$ (in units of $J$). a, Loschmidt zeros for the quench to
$J_{z}=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 $J_{z}=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 $\Phi=\pi/2$ only (see App. D for details).

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) using the parameters $B=0$, $J=J_{z}>0$ and $D/J=2,\,3,\,4$.
The ground states of the post-quench Hamiltonians are within the topologically
trivial large-$D$ phase. The same quenches have been explored in Ref. Hagymási
_et al._ (2019) 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 $\Phi$ 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\simeq
1$ (with $J=J_{z}$ and $B=0$) separating the Haldane phase from the large-$D$
phase Chen  _et al._ (2003), providing 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 finite-size effects. In this
prescription, a thermodynamic line of zeros is approximated by the smooth line
of zeros emerging at the twist angle $\Phi=\pi/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. Hagymási  _et al._ (2019). However, in contrast to
Ref. Hagymási  _et al._ (2019), which considers 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 Heyl
(2018), 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 Chen  _et al._ (2003). 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}\simeq 1.2J$
Chen  _et al._ (2003). 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 Haldane-to-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\leq 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.

Figure 5: Determination of the critical time from the initial energy
fluctuations. Loschmidt zeros for the Heisenberg chain (11) obtained from the
energy fluctuations in the ground state of the model for $J=J_{z}=3B$ and
$D=0$ at the initial time $\tau=0$, while the energy is determined by the
post-quench Hamiltonian with $B=0$ and $D=4$. Here $J=J_{z}=1$ is the unit of
energy and inverse time. The zeros correspond to chains of lengths
$L=5,\ldots,9$, and in the upper (lower) panel we have extracted the zeros
using energy cumulants of orders $n=4,\ldots,14$ ($n=8,\ldots,19$).
Importantly, the zeros converge to their exact positions with increasing
cumulant orders as it can be seen by comparing the panels. The grey line
corresponds to the zeros in panel c1 of Fig. 3, and the estimate of the
critical time is indicated with a red cross.

The Loschmidt moments are generally complex-valued, and it is not obvious how
they can be measured. However, at the initial time, $\tau=0$, the Loschmidt
moments simplify to the moments of the post-quench Hamiltonian with respect to
the initial state as
$\langle\mathcal{\hat{H}}^{n}\rangle_{0}=\left\langle\Psi_{0}\right|\mathcal{\hat{H}}^{n}\left|\Psi_{0}\right\rangle$.
Thus, by repeatedly preparing the system in the state
$\left|\Psi_{0}\right\rangle$ and measuring the energy given by the post-
quench Hamiltonian $\mathcal{\hat{H}}$, 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,\ldots,9$. From these results, we predict the critical time to be around
$t_{c}\simeq 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 Eckardt (2017) or
with spin chains on surfaces Choi  _et al._ (2019), 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 Flindt  _et al._ (2009);
Fricke  _et al._ (2010) and used them for determining the zeros of generating
functions Flindt and Garrahan (2013); Brandner  _et al._ (2017), 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 near-term quantum computers, which
are now becoming available. Such quantum computers allow for the specific
tailoring of any desired Hamiltonian and for time-evolving an initial state
both in real and imaginary time Sun  _et al._ (2021); Lin  _et al._ (2021).
Thus, it will be possible to evaluate a time-evolved state of the form
$\left|\Psi(\tau)\right\rangle=e^{-\tau\mathcal{\hat{H}}}\left|\Psi_{0}\right\rangle$
and subsequently calculate the Loschmidt moments
$\langle\mathcal{\hat{H}}^{n}\rangle_{\tau}=\left\langle\Psi_{0}\right|\mathcal{\hat{H}}^{n}\left|\Psi(\tau)\right\rangle/\langle\Psi_{0}|\Psi(\tau)\rangle$
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.

###### Acknowledgements.

We thank the authors of Ref. Hagymási  _et al._ (2019) for providing us with
their results for the spin-1 Heisenberg chain, which we used to extract the
critical times indicated in Fig. 3. The work was supported by the Academy of
Finland through the Finnish Centre of Excellence in Quantum Technology
(project numbers 312057 and 312299) as well as projects number 308515, 330384,
331094, and 331737. F.B. acknowledges support from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement number 892956. T.O. acknowledges project funding from
Helsinki Institute of Physics.

##  Appendix A Determination of Loschmidt zeros

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
($\tau_{k}\neq\tau_{k^{\prime}}$ for $k\neq k^{\prime}$), and also to denote
by $d_{k}$ the multiplicity of the $k$-th zero. Using this convention, Eq. (7)
becomes

$\kappa_{n}(\tau)=\sum_{k=0}^{\infty}\frac{d_{k}}{(\tau_{k}-\tau)^{n}}\simeq\sum^{m-1}_{k=0}d_{k}\lambda_{k}^{n}\,.$
(12)

On the right hand side, we have introduced $\lambda_{k}=1/(\tau_{k}-\tau)$ 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 $\kappa_{n}(\tau)$ 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. (12) were exact, the normalized cumulants would satisfy a homogeneous
linear difference equation of degree $m$ of the form

$\kappa_{n}=a_{1}\kappa_{n-1}+a_{2}\kappa_{n-2}+\ldots+a_{m}$ (13)

for some coefficients $a_{l}$. Indeed, the general solution of Eq. (13) is
given by the right hand side of Eq. (12), where $\lambda_{k}$ are the
characteristic roots of the characteristic equation associated to Eq. (13),
see Eq. (15) below, and $d_{k}$ can be arbitrary coefficients due to
linearity. This observation is crucial for inverting Eq. (12) and extracting
the zeros from the cumulants as we explain below. We note that Eqs. (12,13)
are exact only if the Loschmidt amplitude $\mathcal{Z}(\tau)$ is a polynomial
with $m$ distinct zeros. In this case, the method provides exactly all the
zeros of $\mathcal{Z}(\tau)$, independently of the cumulant orders used. In
the general case, where $\mathcal{Z}(\tau)$ is an entire function, the method
provides approximate zeros $\tau_{k}^{(n)}$ that depend on the cumulant orders
and converge to the exact zeros for $n\to\infty$. Associated to the
approximate zeros, the method provides also a sequence of approximate
multiplicities $d_{k}^{(n)}$ 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,\dots,m}$
in Eq. (13). This can be done by solving a linear system of $m$ equations,
which requires the knowledge of $2m$ consecutive cumulants ($\kappa_{l}$, with
$n-m\leq l\leq n+m-1$) and takes the form

$\begin{pmatrix}\kappa_{n-1}&\kappa_{n-2}&\dots&\kappa_{n-m+1}&\kappa_{n-m}\\\
\kappa_{n}&\kappa_{n-1}&\dots&\kappa_{n-m+2}&\kappa_{n-m+1}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
\kappa_{n+m-3}&\kappa_{n+m-4}&\dots&\kappa_{n-1}&\kappa_{n-2}\\\
\kappa_{n+m-2}&\kappa_{n+m-3}&\dots&\kappa_{n}&\kappa_{n-1}\end{pmatrix}\begin{pmatrix}a^{(n)}_{1}\\\
a^{(n)}_{2}\\\ \vdots\\\ a^{(n)}_{m-1}\\\
a^{(n)}_{m}\end{pmatrix}=\begin{pmatrix}\kappa_{n}\\\ \kappa_{n+1}\\\
\vdots\\\ \kappa_{n+m-2}\\\ \kappa_{n+m-1}\end{pmatrix}.$ (14)

The square matrix on the left hand side is a Toeplitz matrix. With the
notation $a_{l}^{(n)}$, we emphasize that the coefficients of the linear
difference equation obtained from Eq. (14) depend on the cumulant orders used,
since Eqs. (12,13) are approximations in general.

The second step is to solve the characteristic equation associated to the
linear difference equation (13),

$\lambda^{m}-a^{(n)}_{1}\lambda^{m-1}-a^{(n)}_{2}\lambda^{m-2}-\ldots-a^{(n)}_{m-1}\lambda-a^{(n)}_{m}=0,$
(15)

which is a polynomial equation in $\lambda$. This provides the $m$
characteristic roots $\lambda^{(n)}_{k=0,\dots,m-1}$ of the difference
equation (13). Then, the approximate zeros are obtained from the relation
$\tau^{(n)}_{k}=\tau+1/\lambda^{(n)}_{k}$. In Eq. (14), we have suppressed the
dependence of the cumulants $\kappa_{n}(\tau)$ on the basepoint $\tau$ at
which they are calculated, see Eqs. (7, 12). However, we emphasize that the
approximate zeros $\tau^{(n)}_{k}$ depend in general not only on the cumulant
orders used in Eq. (14), but also on the basepoint $\tau$. As an example, by
solving explicitly Eqs. (14,15) in the case $m=2$, one obtains

$\begin{split}(\tau^{(n)}_{0}-\tau)+(\tau^{(n)}_{1}-\tau)&=\frac{\kappa_{n}\kappa_{n-1}-\kappa_{n-2}\kappa_{n+1}}{\kappa_{n}^{2}-\kappa_{n+1}\kappa_{n-1}},\\\
(\tau_{0}^{(n)}-\tau)(\tau_{1}^{(n)}-\tau)&=\frac{\kappa_{n-1}^{2}-\kappa_{n}\kappa_{n-2}}{\kappa_{n}^{2}-\kappa_{n+1}\kappa_{n-1}},\end{split}$
(16)

where $\tau_{0}^{(n)}$ and $\tau_{1}^{(n)}$ are two sequences that converge
for $n\to\infty$ to the two closest zeros $\tau_{0}$ and $\tau_{1}$
($\neq\tau_{0}^{*}$ in general).

The third and final step is the calculation of the coefficients $d_{k}^{(n)}$
by solving the following linear system of equations,

$\begin{pmatrix}1&1&\dots&1&1\\\
\lambda_{0}&\lambda_{1}&\dots&\lambda_{m-2}&\lambda_{m-1}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
\lambda_{0}^{m-2}&\lambda_{1}^{m-2}&\dots&\lambda_{m-2}^{m-2}&\lambda_{m-1}^{m-2}\\\\[3.00003pt]
\lambda_{0}^{m-1}&\lambda_{1}^{m-1}&\dots&\lambda_{m-2}^{m-1}&\lambda_{m-1}^{m-1}\end{pmatrix}\begin{pmatrix}d^{(n)}_{0}\lambda_{0}^{n}\\\
d^{(n)}_{1}\lambda_{1}^{n}\\\ \vdots\\\
d^{(n)}_{m-2}\lambda_{m-2}^{n}\\\\[3.00003pt]
d^{(n)}_{m-1}\lambda_{m-1}^{n}\end{pmatrix}=\begin{pmatrix}\kappa_{n}\\\
\kappa_{n+1}\\\ \vdots\\\ \kappa_{n+m-2}\\\\[3.00003pt]
\kappa_{n+m-1}\end{pmatrix}\,,$ (17)

where we have dropped the superscript from $\lambda_{k}^{(n)}$ in the above
equation to ease the notation. The matrix on the left-hand side is a
Vandermonde matrix, which is invertible, if all the $\lambda^{(n)}_{k}$ are
distinct. Again, we emphasize that the approximate multiplicities
$d_{k}^{(n)}$ depend on the chosen cumulant orders and on the basepoint,
moreover, they are not exactly integers in general since Eqs. (12,13) are
approximations. The remarkable property of the coefficients $d_{k}^{(n)}$
obtained from Eq. (17) is that they converge to the respective multiplicities
$d_{k}^{(n)}\to d_{k}$ for $n\to\infty$ together with the approximate zeros
$\tau_{k}^{(n)}\to\tau_{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 zeros - numerical 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. (14)) to observe the
convergence of the approximate zeros $\tau_{k}^{(n)}$. 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. (17) converge to the exact multiplicities $d_{k}^{(n)}\to
d_{k}$ concomitantly with the approximate zeros $\tau_{n}^{(n)}\to\tau_{k}$.
Indeed, in the case where the exact zeros $\tau_{k}$ and their multiplicities
$d_{k}$ are known in advance (for instance the Kitaev chain with $V=0$), we
have observed that the distance $|\tau_{k}^{(n)}-\tau_{k}|$ is approximately
linearly proportional to $|d_{k}^{(n)}-d_{k}|$. This is an extremely useful
fact that allows us to automatically select the $\tau_{k}^{(n)}$ which are
good approximations of the exact zeros. This is done by retaining only the
pairs $(\tau_{k}^{(n)},\,d_{k}^{(n)})$ for which $|d_{k}^{(n)}-\ell|<r$, where
$r$ is a fixed threshold and $\ell$ is a chosen integer. This condition
appears to be necessary but not sufficient to guarantee that $\tau_{k}^{(n)}$
is a good approximation of $\tau_{k}$. Indeed, it can occur that the
approximate zero of a pair $(\tau_{k}^{(n)},d_{k}^{(n)})$ is not a good
approximation of any exact zero even if the condition $|d^{(n)}_{k}-\ell|<r$
is satisfied for some integer $\ell$. 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
($\tau_{0}=\tau_{1}^{*}$, $m=2$) has been applied earlier to study critical
phenomena in classical equilibrium problems Deger  _et al._ (2018); Deger and
Flindt (2019, 2020); Deger  _et al._ (2020), 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_{k}^{(n)}$ 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 $\tau$-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 $\tau=0$ and ending at
$\tau=-4+6i$. In this way we can map out zeros in different regions of the
complex plane. The intermediate values of $\tau$ 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 $\delta\tau=0.001$. For each
basepoint $\tau$ we apply Eqs. (14), (15) and (17) using cumulants
$\kappa_{n}(\tau)$ from $n=9$ to $n=22$, and we obtain $m=7$ distinct
approximate zeros $\tau_{k=0,\dots,6}^{(j)}$ and corresponding coefficients
$d_{k=0.,\dots,6}^{(j)}$. Here we change the notation slightly: the
superscript in the pair $(\tau_{k}^{(j)},d_{k}^{(j)})$ 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 $(\tau_{k}^{(j)},d_{k}^{(j)})$
for which $|d_{k}^{(j)}-1|<r=0.01$. We have not found zeros with multiplicity
$\ell>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$ Heyl
_et al._ (2013), 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
Arthur and Vassilvitskii (2007) 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
$(\tau_{k}^{(j)},d_{k}^{(j)})$. Clusters with less then ten pairs are
discarded to eliminate any false positives.

Finally, within each cluster we select the pair $(\tau_{k}^{(j)},d_{k}^{(j)})$
for which $|d_{k}^{(j)}-1|$ takes its smallest value. The approximate zero
$\tau_{k}^{(j)}$ 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 $|\tau_{k}^{(j)}-\tau_{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\approx 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. 2a-c. 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 error. 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 $\tau$-plane is performed
by using the standard Krylov subspace method Park and Light (1986); Paeckel
_et al._ (2019). A time step $\delta\tau$ in the evolution is performed by
first computing an orthonormal basis
$B=\\{\left|v_{0}\right\rangle=\left|\Psi(\tau)\right\rangle,\left|v_{1}\right\rangle,\dots,\left|v_{N_{\rm
vec}}\right\rangle\\}$ of the Krylov subspace

$\mathcal{K}_{N_{\rm
vec}}=\mathrm{Span}\\{\mathcal{\hat{H}}^{n}\left|\Psi(\tau)\right\rangle|\,\,n=0,\dots,N_{\rm
vec}\\}\,.$ (18)

We use the QR decomposition to compute the orthonormal basis of
$\mathcal{K}_{N_{\rm vec}}$ to ensure that there is no loss of orthogonality
as in the standard Lanczos algorithm Paeckel  _et al._ (2019). Then, the
approximate time-evolved state is obtained as

$\left|\Psi(\tau+\delta\tau)\right\rangle\simeq e^{-\delta\tau\hat{H}_{\rm
eff}}\left|\Psi(\tau)\right\rangle\,,$ (19)

where $\hat{H}_{\rm eff}$ is the effective Hamiltonian, that is an operator
which acts on the Krylov subspace and whose matrix elements are given by
$\left\langle v_{i}\right|\mathcal{\hat{H}}\left|v_{j}\right\rangle$ with
$\left|v_{i}\right\rangle\in B$. It is represented by a square matrix of
dimension $N_{\rm vec}+1$ whose exponential is easy to evaluate since we take
$N_{\rm vec}=8$ in our case. As suggested in Section 5 of Ref. Paeckel  _et
al._ (2019), the effective Hamiltonian is forcibly set to be a tridiagonal
matrix to improve stability.

In order to estimate the accuracy of Eq. (19) we perform the time evolution on
two Kyrlov subspaces $\mathcal{K}_{N^{\prime}_{\rm vec}}$ and
$\mathcal{K}_{N^{\prime\prime}_{\rm vec}}$ with $N^{\prime}_{\rm
vec}+1=N^{\prime\prime}_{\rm vec}\leq N_{\rm vec}$ and compute the distance
between the approximate evolved states
$d=\|\left|\Psi^{\prime}(\tau+\delta\tau)\right\rangle-\left|\Psi^{\prime\prime}(\tau+\delta\tau)\right\rangle\|_{2}$.
If $d<10^{-10}$ the state
$\left|\Psi^{\prime\prime}(\tau+\delta\tau)\right\rangle$ is stored and used
to evaluate the moments of the Loschmidt amplitude according to
$\langle\mathcal{\hat{H}}^{n}\rangle_{\tau+\delta\tau}=\left\langle\Psi_{0}\right|\mathcal{\hat{H}}^{n}\left|\Psi^{\prime\prime}(\tau+\delta\tau)\right\rangle$.
On the other hand, if for a given $\delta\tau$ the condition is not satisfied
even for $N^{\prime\prime}_{\rm vec}=N_{\rm vec}$, the Krylov subspace
$\mathcal{K}_{N_{\rm vec}}$ in Eq. (18) is recomputed using as the seed state
the last stored vector
$\left|\Psi^{\prime\prime}(\tau+\delta\tau)\right\rangle$. The whole process
is repeated up to the desired final $\tau$.

An important advantage of the Krylov time evolution algorithm described above
is that intermediate values of $\tau$ 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 (18). In order to compute the moments and the cumulants, a
larger Krylov subspace ($N_{\rm vec}=22$ in our case) has to be calculated
only at $\tau=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 Fáth and Sólyom (1993). Here we discuss how
these boundary conditions can also help to gauge finite-size effects in the
Loschmidt zeros. For the Haldane chain, Eq. (11), the spectrum exhibits a
$1/L$ dependence on the boundary conditions Fáth and Sólyom (1993). 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 Torre  _et al._
. On the other hand, in finite systems, the many-body energies $E_{i}(\Phi)$
depend on the twist angle $\Phi$, in particular the highly excited states. The
Loschmidt zeros inherit this property – a single zero, say $\tau_{k}$, gives
rise to a multiplet of zeros $\tau_{k}(\Phi)$ corresponding to different
values of $\Phi$ 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}(\Phi)$ are even functions of $\Phi$,
having their extreme values at $\Phi=0$ or $\Phi=\pi$. As seen in Fig. 3a1,
the boundary points of the multiplets of zeros exactly correspond to $\Phi=0$
or $\Phi=\pi$. In the thermodynamic limit, these multiplets converge to a
single point, $\tau_{k}$, for all angles $\Phi$.

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 $\Phi=0$ and $\pi$. The
curve of zeros corresponding to a general value of $\Phi$ exhibits size-
dependent oscillations around the mean value as seen in Fig. 3a1 by the zig-
zag line corresponding to $\Phi=0$. However, these finite-size oscillations
are suppressed for zeros corresponding to the twist angle $\Phi=\pi/2$, which
form a smooth curve. This observation suggests that the twist angle
$\Phi=\pi/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. Hagymási  _et al._ (2019), which are also indicated in Fig. 3.
A similar pattern is observed for the quench from the Haldane to the Néel
phase in Fig. 4, providing further support for approximating the thermodynamic
lines of zeros using the zeros corresponding to the boundary condition
$\Phi=\pi/2$.

## References

  * Chandler (1987) D. Chandler,  _Introduction to Modern Statistical Mechanics_ (Oxford University Press, New York, 1987).
  * Goldenfeld (2018) N. Goldenfeld,  _Lectures on Phase Transitions and the Renormalization Group_ (CRC Press, Boca Raton, 2018).
  * Heyl  _et al._ (2013) M. Heyl, A. Polkovnikov, and S. Kehrein, “Dynamical Quantum Phase Transitions in the Transverse-Field Ising Model,” Phys. Rev. Lett. 110, 135704 (2013).
  * Karrasch and Schuricht (2013) C. Karrasch and D. Schuricht, “Dynamical phase transitions after quenches in nonintegrable models,” Phys. Rev. B 87, 195104 (2013).
  * Andraschko and Sirker (2014) F. Andraschko and J. Sirker, “Dynamical quantum phase transitions and the loschmidt echo: A transfer matrix approach,” Phys. Rev. B 89, 125120 (2014).
  * Vajna and Dóra (2014) S. Vajna and B. Dóra, “Disentangling dynamical phase transitions from equilibrium phase transitions,” Phys. Rev. B 89, 161105 (2014).
  * Vajna and Dóra (2015) S. Vajna and B. Dóra, “Topological classification of dynamical phase transitions,” Phys. Rev. B 91, 155127 (2015).
  * Budich and Heyl (2016) J. C. Budich and M. Heyl, “Dynamical topological order parameters far from equilibrium,” Phys. Rev. B 93, 085416 (2016).
  * Zvyagin (2016) A. A. Zvyagin, “Dynamical quantum phase transitions (Review Article),” J. Low Temp. Phys. 42, 971 (2016).
  * Halimeh and Zauner-Stauber (2017) J. C. Halimeh and V. Zauner-Stauber, “Dynamical phase diagram of quantum spin chains with long-range interactions,” Phys. Rev. B 96, 134427 (2017).
  * Heyl (2018) M. Heyl, “Dynamical quantum phase transitions: a review,” Rep. Prog. Phys. 81, 054001 (2018).
  * Jurcevic  _et al._ (2017) P. Jurcevic, H. Shen, P. Hauke, C. Maier, T. Brydges, C. Hempel, B. P. Lanyon, M. Heyl, R. Blatt, and C. F. Roos, “Direct Observation of Dynamical Quantum Phase Transitions in an Interacting Many-Body System,” Phys. Rev. Lett. 119, 080501 (2017).
  * Zhang  _et al._ (2017) J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, “Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator,” Nature 551, 601 (2017).
  * Fläschner  _et al._ (2018) N. Fläschner, D. Vogel, M. Tarnowski, B. S. Rem, D.-S. Lühmann, M. Heyl, J. C. Budich, L. Mathey, K. Sengstock, and C. Weitenberg, “Observation of dynamical vortices after quenches in a system with topology,” Nat. Phys. 14, 265 (2018).
  * Guo  _et al._ (2019) X.-Y. Guo, C. Yang, Y. Zeng, Y. Peng, H.-K. Li, H. Deng, Y.-R. Jin, S. Chen, D. Zheng, and H. Fan, “Observation of a Dynamical Quantum Phase Transition by a Superconducting Qubit Simulation,” Phys. Rev. Applied 11, 044080 (2019).
  * Tian  _et al._ (2019) T. Tian, Y. Ke, L. Zhang, S. Lin, Z. Shi, P. Huang, C. Lee, and J. Du, “Observation of dynamical phase transitions in a topological nanomechanical system,” Phys. Rev. B 100, 024310 (2019).
  * Wang  _et al._ (2019) K. Wang, X. Qiu, L. Xiao, X. Zhan, Z. Bian, W. Yi, and P. Xue, “Simulating Dynamic Quantum Phase Transitions in Photonic Quantum Walks,” Phys. Rev. Lett. 122, 020501 (2019).
  * Tian  _et al._ (2020) T. Tian, H.-X. Yang, L.-Y. Qiu, H.-Y. Liang, Y.-B. Yang, Y. Xu, and L.-M. Duan, “Observation of Dynamical Quantum Phase Transitions with Correspondence in an Excited State Phase Diagram,” Phys. Rev. Lett. 124, 043001 (2020).
  * Xu  _et al._ (2020) K. Xu, Z.-H. Sun, W. Liu, Y.-R. Zhang, H. Li, H: Dong, W. Ren, P. Zhang, F. Nori, D. Zheng, H. Fan, and H. Wang, “Probing dynamical phase transitions with a superconducting quantum simulator,” Science Adv. 6, eaba4935 (2020).
  * Pozsgay (2013) B. Pozsgay, “The dynamical free energy and the Loschmidt echo for a class of quantum quenches in the Heisenberg spin chain,” J. Stat. Mech. 2013, P10028 (2013).
  * Kriel  _et al._ (2014) J. N. Kriel, C. Karrasch, and S. Kehrein, “Dynamical quantum phase transitions in the axial next-nearest-neighbor Ising chain,” Phys. Rev. B 90, 125106 (2014).
  * Sharma  _et al._ (2015) S. Sharma, S. Suzuki, and A. Dutta, “Quenches and dynamical phase transitions in a nonintegrable quantum Ising model,” Phys. Rev. B 92, 104306 (2015).
  * Zauner-Stauber and Halimeh (2017) V. Zauner-Stauber and J. C. Halimeh, “Probing the anomalous dynamical phase in long-range quantum spin chains through Fisher-zero lines,” Phys. Rev. E 96, 062118 (2017).
  * Homrighausen  _et al._ (2017) I. Homrighausen, N. O. Abeling, V. Zauner-Stauber, and J. C. Halimeh, “Anomalous dynamical phase in quantum spin chains with long-range interactions,” Phys. Rev. B 96, 104436 (2017).
  * Heyl  _et al._ (2018) M. Heyl, F. Pollmann, and B. Dóra, “Detecting Equilibrium and Dynamical Quantum Phase Transitions in Ising Chains via Out-of-Time-Ordered Correlators,” Phys. Rev. Lett. 121, 016801 (2018).
  * Kennes  _et al._ (2018) D. M. Kennes, D. Schuricht, and C. Karrasch, “Controlling dynamical quantum phase transitions,” Phys. Rev. B 97, 184302 (2018).
  * Hagymási  _et al._ (2019) I. Hagymási, C. Hubig, Ö. Legeza, and U. Schollwöck, “Dynamical Topological Quantum Phase Transitions in Nonintegrable Models,” Phys. Rev. Lett. 122, 250601 (2019).
  * Lacki and Heyl (2019) M. Lacki and M. Heyl, “Dynamical quantum phase transitions in collapse and revival oscillations of a quenched superfluid,” Phys. Rev. B 99, 121107 (2019).
  * Huang  _et al._ (2019) Y.-P. Huang, D. Banerjee, and M. Heyl, “Dynamical Quantum Phase Transitions in U(1) Quantum Link Models,” Phys. Rev. Lett. 122, 250401 (2019).
  * Halimeh  _et al._ (2020) J. C. Halimeh, M. Van Damme, V. Zauner-Stauber, and L. Vanderstraeten, “Quasiparticle origin of dynamical quantum phase transitions,” Phys. Rev. Research 2, 033111 (2020).
  * Žunkovič  _et al._ (2018) B. Žunkovič, M. Heyl, M. Knap, and A. Silva, “Dynamical Quantum Phase Transitions in Spin Chains with Long-Range Interactions: Merging Different Concepts of Nonequilibrium Criticality,” Phys. Rev. Lett. 120, 130601 (2018).
  * Weidinger  _et al._ (2017) S. A. Weidinger, M. Heyl, A. Silva, and M. Knap, “Dynamical quantum phase transitions in systems with continuous symmetry breaking,” Phys. Rev. B 96, 134313 (2017).
  * Feldmeier  _et al._ (2019) J. Feldmeier, F. Pollmann, and M. Knap, “Emergent Glassy Dynamics in a Quantum Dimer Model,” Phys. Rev. Lett. 123, 040601 (2019).
  * Heyl (2015) M. Heyl, “Scaling and Universality at Dynamical Quantum Phase Transitions,” Phys. Rev. Lett. 115, 140602 (2015).
  * Schmitt and Kehrein (2015) M. Schmitt and S. Kehrein, “Dynamical quantum phase transitions in the Kitaev honeycomb model,” Phys. Rev. B 92, 075114 (2015).
  * Fogarty  _et al._ (2017) T. Fogarty, A. Usui, T. Busch, A. Silva, and J. Goold, “Dynamical phase transitions and temporal orthogonality in one-dimensional hard-core bosons: from the continuum to the lattice,” New J. Phys. 19, 113018 (2017).
  * Bhattacharya and Dutta (2017) U. Bhattacharya and A. Dutta, “Emergent topology and dynamical quantum phase transitions in two-dimensional closed quantum systems,” Phys. Rev. B 96, 014302 (2017).
  * Sedlmayr  _et al._ (2018) N. Sedlmayr, P. Jaeger, M. Maiti, and J. Sirker, “Bulk-boundary correspondence for dynamical phase transitions in one-dimensional topological insulators and superconductors,” Phys. Rev. B 97, 064304 (2018).
  * Jafari  _et al._ (2019) R. Jafari, H. Johannesson, A. Langari, and M. A. Martin-Delgado, “Quench dynamics and zero-energy modes: The case of the Creutz model,” Phys. Rev. B 99, 054302 (2019).
  * Defenu  _et al._ (2019) N. Defenu, T. Enss, and J. C. Halimeh, “Dynamical criticality and domain-wall coupling in long-range Hamiltonians,” Phys. Rev. B 100, 014434 (2019).
  * Najafi  _et al._ (2019) K. Najafi, M. A. Rajabpour, and J. Viti, “Return amplitude after a quantum quench in the XY chain,” J. Stat. Mech. 2019, 083102 (2019).
  * Gulácsi  _et al._ (2020) B. Gulácsi, M. Heyl, and B. Dóra, “Geometrical quench and dynamical quantum phase transition in the $\alpha-{T}_{3}$ lattice,” Phys. Rev. B 101, 205135 (2020).
  * Zamani  _et al._ (2020) S. Zamani, R. Jafari, and A. Langari, “Floquet dynamical quantum phase transition in the extended xy model: Nonadiabatic to adiabatic topological transition,” Phys. Rev. B 102, 144306 (2020).
  * Deger  _et al._ (2018) A. Deger, K. Brandner, and C. Flindt, “Lee-Yang zeros and large-deviation statistics of a molecular zipper,” Phys. Rev. E 97, 012115 (2018).
  * Deger and Flindt (2019) A. Deger and C. Flindt, “Determination of universal critical exponents using Lee-Yang theory,” Phys. Rev. Research 1, 023004 (2019).
  * Deger and Flindt (2020) A. Deger and C. Flindt, “Lee-Yang theory of the Curie-Weiss model and its rare fluctuations,” Phys. Rev. Research 2, 033009 (2020).
  * Deger  _et al._ (2020) A. Deger, F. Brange, and C. Flindt, “Lee-Yang theory, high cumulants, and large-deviation statistics of the magnetization in the Ising model,” Phys. Rev. B 102, 174418 (2020).
  * Yang and Lee (1952) C. N. Yang and T. D. Lee, “Statistical Theory of Equations of State and Phase Transitions. I. Theory of Condensation,” Phys. Rev. 87, 404 (1952).
  * Lee and Yang (1952) T. D. Lee and C. N. Yang, “Statistical Theory of Equations of State and Phase Transitions. II. Lattice Gas and Ising Model,” Phys. Rev. 87, 410 (1952).
  * Blythe and Evans (2003) R. A. Blythe and M. R. Evans, “The Lee-Yang theory of equilibrium and nonequilibrium phase transitions,” Braz. J. Phys. 33, 464 (2003).
  * Bena  _et al._ (2005) I. Bena, M. Droz, and A. Lipowski, “Statistical mechanics of equilibrium and nonequilibrium phase transitions: the Yang-Lee formalism,” Int. J. Mod. Phys. B 19, 4269 (2005).
  * Arfken  _et al._ (2012) G. B. Arfken, H. J. Weber, and F. E. Harris,  _Mathematical Methods for Physicists: A Comprehensive Guide_ (Academic Press, New York, 2012).
  * (53) M. E. Fisher, in  _Boulder Lectures in Theoretical Physics_ , Vol. 7 (University of Colorado, Boulder).
  * Park and Light (1986) T. J. Park and J. C. Light, “Unitary quantum time evolution by iterative Lanczos reduction,” J. Chem. Phys. 85, 5870 (1986).
  * Paeckel  _et al._ (2019) S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, “Time-evolution methods for matrix-product states,” Ann. Phys. 411, 167998 (2019).
  * Franchini (2016) F. Franchini,  _An introduction to integrable techniques for one-dimensional quantum systems_ (Springer, 2016).
  * Affleck  _et al._ (1987) I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, “Rigorous results on valence-bond ground states in antiferromagnets,” Phys. Rev. Lett. 59, 799 (1987).
  * Chen  _et al._ (2003) W. Chen, K. Hida, and B. C. Sanctuary, “Ground-state phase diagram of $S=1$ $\mathrm{XXZ}$ chains with uniaxial single-ion-type anisotropy,” Phys. Rev. B 67, 104401 (2003).
  * Eckardt (2017) A. Eckardt, “Colloquium: Atomic quantum gases in periodically driven optical lattices,” Rev. Mod. Phys. 89, 011004 (2017).
  * Choi  _et al._ (2019) D.-J. Choi, N. Lorente, J. Wiebe, K. von Bergmann, A. F. Otte, and A. J. Heinrich, “Colloquium: Atomic spin chains on surfaces,” Rev. Mod. Phys. 91, 041001 (2019).
  * Flindt  _et al._ (2009) C. Flindt, C. Fricke, F. Hohls, T. Novotný, K. Netočný, T. Brandes, and R. J. Haug, “Universal oscillations in counting statistics,” Proc. Natl. Acad. Sci. USA 106, 10116 (2009).
  * Fricke  _et al._ (2010) C. Fricke, F. Hohls, C. Flindt, and R. J. Haug, “High cumulants in the counting statistics measured for a quantum dot,” Physica E 42, 848 (2010).
  * Flindt and Garrahan (2013) C. Flindt and J. P. Garrahan, “Trajectory Phase Transitions, Lee-Yang Zeros, and High-Order Cumulants in Full Counting Statistics,” Phys. Rev. Lett. 110, 050601 (2013).
  * Brandner  _et al._ (2017) K. Brandner, V. F. Maisi, J. P. Pekola, J. P. Garrahan, and C. Flindt, “Experimental Determination of Dynamical Lee-Yang Zeros,” Phys. Rev. Lett. 118, 180601 (2017).
  * Sun  _et al._ (2021) S.-N. Sun, M. Motta, R. N. Tazhigulov, A. T.K. Tan, G. K.-L. Chan, and A. J. Minnich, “Quantum Computation of Finite-Temperature Static and Dynamical Properties of Spin Systems Using Quantum Imaginary Time Evolution,” Phys. Rev. X Quantum 2, 010317 (2021).
  * Lin  _et al._ (2021) S.-H. Lin, R. Dilip, A. G. Green, A. Smith, and F. Pollmann, “Real- and Imaginary-Time Evolution with Compressed Quantum Circuits,” Phys. Rev. X Quantum 2, 010342 (2021).
  * Arthur and Vassilvitskii (2007) D. Arthur and S. Vassilvitskii, “$k$-means++: the advantages of careful seeding,” in  _Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms_ (Society for Industrial and Applied Mathematics Philadelphia, PA, USA, 2007) p. 1027.
  * Fáth and Sólyom (1993) G. Fáth and J. Sólyom, “Isotropic spin-1 chain with twisted boundary condition,” Phys. Rev. B 47, 872 (1993).
  * (69) G. Torre, V. Marić, D. Kuić, F. Franchini, and S. M. Giampaolo, “An odd thermodynamic limit for the Loschmidt echo,” arXiv:2105.06483 .