# Detecting quantum phase transitions in the quasi-stationary regime of Ising
chains

Ceren B. Dağ<EMAIL_ADDRESS>ITAMP, Harvard-Smithsonian Center for
Astrophysics, Cambridge, Massachusetts, 02138, USA  Department of Physics,
Harvard University, 17 Oxford Street Cambridge, MA 02138, USA  Department of
Physics, University of Michigan, Ann Arbor, Michigan 48109, USA  Philipp
Uhrich  INO-CNR BEC Center and Department of Physics, University of Trento,
Via Sommarive 14, I-38123 Trento, Italy  Yidan Wang  Department of Physics,
Harvard University, 17 Oxford Street Cambridge, MA 02138, USA  Ian P.
McCulloch  School of Mathematics and Physics, University of Queensland, St.
Lucia, Queensland 4072, Australia  Jad C. Halimeh  INO-CNR BEC Center and
Department of Physics, University of Trento, Via Sommarive 14, I-38123 Trento,
Italy

###### Abstract

Recently, single-site observables have been shown to be useful for the
detection of dynamical criticality due to an emergence of a universal
critically-prethermal temporal regime in the magnetization [arXiv:2105.05986].
Here, we explore the potential of single-site observables as probes of quantum
phase transitions in integrable and nonintegrable transverse-field Ising
chains (TFIC). We analytically prove the requirement of zero modes for a
quasi-stationary temporal regime to emerge at a bulk probe site, and show how
this regime gives rise to a non-analytic behavior in the dynamical order
profiles. Our $t$-DMRG calculations verify the results of the quench mean-
field theory for near-integrable TFIC both with finite-size and finite-time
scaling analyses. We find that both finite-size and finite-time analyses
suggest a dynamical critical point for a strongly nonintegrable and locally
connected TFIC. We finally demonstrate the presence of a quasi-stationary
regime in the power-law interacting TFIC, and extract local dynamical order
profiles for TFIC in the long-range Ising universality class with algebraic
light cones.

###### Contents

  1. I Introduction
  2. II The quasi-stationary regime in the short-range TFIC
  3. III The integrable TFIC
    1. III.1 Quench dynamics of single-site magnetization via the cluster theorem
    2. III.2 Requirement of zero modes to observe a quasi-stationary regime
    3. III.3 Derivation of the series expression for the edge magnetization
    4. III.4 Detection of QCP
  4. IV The nonintegrable TFIC with next nearest neighbor coupling
    1. IV.1 Quench Mean-field theory (qMFT) analysis
    2. IV.2 Detection of QCP in the near integrable TFIC through qMFT analysis
    3. IV.3 $t$-DMRG calculations
  5. V The quasi-stationary regime in long-range interacting nonintegrable TFIC
  6. VI Conclusions and Outlook
  7. A Independency of the results from temporal cutoffs
  8. B Initial state independence of the results
  9. C Energy gap analysis
  10. D Changing the boundary conditions
  11. E Benchmarking the qMFT method

##  I Introduction

Phase transitions of matter are among the most prevalent physical processes in
nature, during which certain properties of matter can change discontinuously
Landau (1937); Cardy  _et al._ (1996). The notions of universality and scaling
allow a unified description of microscopically different systems by means of
universality classes, whereby the systems in the same class exhibit identical
critical behavior independently of their microscopic details Kadanoff (1993);
Zinn-Justin (2002); Sachdev (2001). Far from equilibrium, universality and
scaling are well-understood in classical systems Hohenberg and Halperin
(1977), unlike in quantum many-body systems, despite of the exploration of
various analogous concepts, e.g., dynamical quantum phase transitions and
criticality Polkovnikov  _et al._ (2011); Mori  _et al._ (2018); Heyl (2018);
Calabrese and Cardy (2006); Calabrese  _et al._ (2011); Kollath  _et al._
(2007); Calabrese  _et al._ (2012); Tsuji  _et al._ (2013); Heyl  _et al._
(2013); Sciolla and Biroli (2013); Chiocchetta  _et al._ (2015); Halimeh and
Zauner-Stauber (2017); Karl  _et al._ (2017); Titum  _et al._ (2019); Gurarie
(2019); Uhrich  _et al._ (2020); Wu (2020a, b); Dağ and Sun (2021); Halimeh
_et al._ (2021); Haldar  _et al._ (2021); Halimeh  _et al._ (2019); Trapin
_et al._ (2021); Heyl  _et al._ (2018); Dağ  _et al._ (2019); Gómez-Ruiz  _et
al._ (2018); Dağ  _et al._ (2020); Sun  _et al._ (2020); Bandyopadhyay  _et
al._ (2021a, b), and the experimental observation of out-of-equilibrium
critical behavior in various quantum synthetic matter (QSM) setups Nicklas
_et al._ (2015); Jurcevic  _et al._ (2017); Zhang  _et al._ (2017); Hauke  _et
al._ (2014); Fläschner  _et al._ (2018); Yang  _et al._ (2019); Muniz  _et
al._ (2020); Nie  _et al._ (2020); Xu  _et al._ (2020); Chen  _et al._ (2020).

With the advent of QSM experiments enjoying a high level of precision and
control Bloch  _et al._ (2008), including single-site addressing techniques
Bakr  _et al._ (2009), the prospect of probing equilibrium and dynamical
quantum phase transitions following quantum quenches has become realistic and
appealing. Although it may sound counterintuitive, extracting equilibrium
criticality through quench dynamics can be an experimentally more viable and
simpler scheme than the challenging procedure of cooling a quantum many-body
system to its ground state. Instead, the system is initialized in an easily
accessible product state, and subsequently quenched through a control
parameter. This method has enabled the extraction of quantum critical points
and universal scaling laws in various models Sciolla and Biroli (2013);
Nicklas  _et al._ (2015); Chiocchetta  _et al._ (2015); Halimeh and Zauner-
Stauber (2017); Karl  _et al._ (2017); Dağ  _et al._ (2018); Heyl  _et al._
(2018); Yang  _et al._ (2019); Titum  _et al._ (2019); Uhrich  _et al._
(2020); Dağ and Sun (2021); Haldar  _et al._ (2021).

In this vein of probing quantum phase transitions and criticality through
quench dynamics, some of us proposed in Ref. Dağ  _et al._ (2021) to use
single-site observables close to an edge in open-boundary and short-range
many-body models. When boundaries are introduced to a suddenly quenched $1$D
transverse-field Ising chain (TFIC) an onset of a quasi-stationary (q.s.)
temporal regime appears Iglói and Rieger (2011), thus interrupting an
otherwise exponential decay to zero in time. Specifically, Ref. Dağ  _et al._
(2021) analytically and numerically revealed a universal and critical
prethermal regime in the vicinity of the quantum critical point (QCP) with
single-site observables, regardless of probe location, initial states and weak
perturbations to the TFIC. This was shown by means of exact analytical and
numerical calculations where the latter is based on a Pfaffian formalism as
well as quench mean-field theory (qMFT). Ref. Dağ  _et al._ (2021) also
derived a scaling function of the dynamical order parameter in the vicinity of
the QCP and showed that the long-lived prethermal regime makes the scaling
function nonlinear. The presence of a nonlinear dynamical scaling function was
found to be independent of the probe site and initial state, and to be robust
to weak integrability breaking.

In this long paper, our theme is the detection of QCP in integrable and
nonintegrable TFIC through single-site observables near the edge of the chain
$r\ll N/2$. We first determine the QCP in the integrable TFIC by locating a
nonanalyticity in the q.s. value of single-site observables. We find
consistently for all sites sufficiently near the edge that the location of
this nonanalyticity coincides with the equilibrium QCP. We also derive an
analytical expression for the edge magnetization in the integrable TFIC
quenched from a polarized state, which is used to characterize the critically
prethermal regime in Ref. Dağ  _et al._ (2021).

We detail the aforementioned Pfaffian formalism behind the exact results on
single-site observables in open-boundary TFICs, and how MFT is implemented in
quench dynamics, e.g., qMFT method. An analytical form of the latter was first
employed in Ref. Titum  _et al._ (2019) for periodic TFIC. Since we lack space
translational symmetry in open-boundary chains, we modify the qMFT method with
the help of the cluster theorem to make it suitable for numerics. By
benchmarking the qMFT method against exact diagonalization (ED) and time-
dependent density matrix renormalization group ($t$-DMRG) results, we show
that the MFT is reliable within times accessible to $t$-DMRG. Then we use our
qMFT method in TFIC with weak interactions, and detect a non-analyticity in
the dynamical order profile slightly shifted from the location of the QCP. By
following Ref. Titum  _et al._ (2019), we call this non-analytic point a
dynamical critical point (DCP).

Furthermore, we perform extensive $t$-DMRG calculations ranging from near-
integrable to strongly nonintegrable and locally connected TFIC, as well as
long-range power-law interacting TFIC. We provide finite-size and finite-time
analyses on the $t$-DMRG data for the near-integrable model up to system sizes
$N=120$, showing that they support the MFT result on the location of the DCP.
Remarkably, evidence of a DCP in the thermodynamic limit, namely a crossing
between the order profiles at different system sizes up to $N=48$, remains in
the $t$-DMRG data for a strongly nonintegrable and locally connected TFIC,
albeit shifted significantly from its corresponding QCP. To our knowledge,
this is the first evidence of a quench measurement of magnetization that
suggests a DCP for the strongly nonintegrable TFIC with ferromagnetic next
nearest neighbor (n.n.n.) couplings $\Delta/J=1$, where $J$ is the fixed
energy scaling set in the calculation. Previously, a dynamical order profile
was generated in Ref. Heyl  _et al._ (2018) with out-of-time-order correlators
at $\Delta/J=0.5$ with no clear evidence of a crossing point between the order
profiles at different system sizes. Ref. Titum  _et al._ (2019) predicted a
DCP, by using two-point correlators, up to $\Delta/J=0.3$ for chains at
$N=25$, and Ref. Haldar  _et al._ (2021) has recently predicted a QCP with the
same correlator, but for antiferromagnetic n.n.n couplings. Additionally, for
this strongly nonintegrable TFIC, we observe that the finite-time analysis of
our $t$-DMRG data underestimates the location of the DCP, thereby shifting it
further away from the QCP. This reveals a breakdown of the validity regime of
the finite-time analysis, and we elaborate on a possible reason of this
observation.

We show the presence of a q.s. temporal regime in TFIC with power-law decaying
interactions, whose equilibrium criticality properties belong either to short-
or long-range Ising universality class with, respectively, $\alpha\geq 3$ and
$\alpha<3$, where $\alpha$ is the power-law exponent of the interaction
strength. We reveal that single-site observables still reflect a difference in
the quench dynamics between the sites near the edges and in the middle when
the model is long-range interacting $\alpha<3$, up until $\alpha=2$ where
there is no longer an onset of a q.s. regime. Upon closer study of the model
with $\alpha=2.5$, which is long-range with algebraic light cones Foss-Feig
_et al._ (2015), we reveal that the local order profiles at different $r$ are
consistent with each other and suggest a crossover.

Because the q.s. regime is the basis of the criticality detection in our
paper, we analytically show why the zero modes are vital in the origin of this
temporal regime. Besides the requirement of zero modes, the origin of the q.s.
regime can also be explained with translational symmetry breaking and making
an asymmetric measurement with respect to the symmetry center of the chain,
e.g., the middle of the chain. This differentiates the two chain edges as
either near or far from the perspective of quasi-particles in the integrable
models, or of the wavefronts in operator spreading in the nonintegrable
models. Hence, the locality of the underlying model, e.g., the linear or
sublinear light cones in the operator spread, rather than the integrability
seems to be a reason for the presence of the q.s. regime. In fact, we
explicitly show instances of q.s. behaviour appearing in strongly
nonintegrable many-body models and its direct application in probing QCP.
Consistently, the q.s. regime does not occur when the model becomes long-range
with logarithmic light cones Hauke and Tagliacozzo (2013). Furthermore, we
find the presence of a q.s. regime even when the hard edges are smoothed out.
This implies that our theoretical predictions could in principle be tested
even in imperfect conditions where hard boundaries might not be present, e.g.,
in trapped atomic and ionic gases.

Although we focus on the dynamics of single-site observables, since we aim to
dynamically probe ferromagnetic order in our paper, we emphasize that our
method could be easily adjusted to probe other quantum orders, e.g., anti-
ferromagnetic order, by increasing the unit cell of the underlying lattice.
Therefore, our work and Ref. Dağ  _et al._ (2021) are the first examples of
spatially minimal measurements that could probe quantum phase transitions and
nonequilibrium criticality.

This paper consists of the following sections. In Sec. II, we describe the
locally-connected TFIC with open boundaries and introduce its q.s. temporal
regime. Then we focus on the detection of QCP in the nearest-neighbor (n.n.)
integrable limit of TFIC where we also explain the formalism behind our
calculations, express the cluster theorem and prove the requirement of zero
modes to observe a q.s. regime in Sec. III. Later in Sec. IV, we detail our
qMFT method and its results, and systematically discuss the $t$-DMRG results
on our locally-connected nonintegrable models. In Sec. V, we focus on the
quench dynamics of single-site observables, and local order profiles in the
power-law interacting TFIC. We conclude in Sec. VI with an outlook.

##  II The quasi-stationary regime in the short-range TFIC

The TFIC with n.n.n. interaction strength $\Delta$ reads

$\displaystyle
H=-J\sum_{1}^{N-1}\sigma_{r}^{z}\sigma_{r+1}^{z}-\Delta\sum_{1}^{N-2}\sigma_{r}^{z}\sigma_{r+2}^{z}+h\sum_{1}^{N}\sigma_{r}^{x},$
(1)

where $\sigma_{r}^{x,z}$ are the Pauli spin matrices on site $r$, $h$ is the
transverse-field strength, $N$ is the length of the open-boundary chain, and
we fix the n.n. coupling strength $J=1$ as the energy scale. The TFIC is a
paradigmatic model of quantum phase transitions, hosting a critical point at
$h_{c}=1$ when $\Delta=0$ which separates a spontaneously broken symmetry
(ferromagnetic) phase from a symmetry preserving (paramagnetic) phase. When
$\Delta>0$, the low-energy properties remain the same, but the QCP shifts to
favor order $h_{c}>1$.

(a)

(b)

Figure 1: (a) The local magnetization $C_{r}(t)$ in the integrable TFIC with
open-boundary (solid-blue) and periodic-boundary (red-dotted and yellow-
squares) conditions given in Eq. (1) after a quench in the transverse-field
strength from $h_{i}=0$ to $h=0.8$. See Secs. III.1 and III.2 for the
numerical methods used in yellow-squares and solid-blue. The red-dotted line
is the analytical equation for the periodic chain derived in Ref. Dağ and Sun
(2021). (b) The local magnetization $C_{r}(t)$ in open-boundary nonintegrable
TFIC ($\Delta>0$) after a quench in the transverse-field strength from
$h_{i}=0$ and calculated with $t-$DMRG.

Let us consider as initial state the ground state $\Ket{\psi_{0}}$ of $H$ at
initial value $h_{i}$ of the transverse-field strength, and then we quench the
latter to a value $h$. While in the case of periodic-boundary conditions the
single-site nonequilibrium response
$C_{r}(t)=\Bra{\psi_{0}}\sigma^{z}_{r}(t)\Ket{\psi_{0}}$ decays exponentially
to zero Calabrese  _et al._ (2012); Dağ and Sun (2021) (see Fig. 1a), open-
boundary conditions give rise to a nonzero-valued equilibrium regime after a
power-law decay when $h_{i}<h\leq h_{c}=1$. This is referred to as the q.s.
regime Iglói and Rieger (2011). Figure 1a shows the q.s. regime of $C_{3}(t)$
in the integrable TFIC quenched from $h_{i}=0$ to $h=0.8$ and at a system size
$N=1500$. While one could argue that the origin of the q.s. regime is because
of the interference between quasi-particles due to the asymmetric location of
the probe site with respect to the symmetry center of the chain that gives
rise to reflection from the closest edge, we show in the next section the
necessity of zero modes to observe a q.s. regime.

We break integrability in Eq. (1) by taking $\Delta>0$ in Fig. 1b for
different physical parameters (see legend). In all cases, our $t-$DMRG
calculations in open-boundary chains (see Sec. IV.3 for the method) show the
emergence of a q.s. regime in relatively early times as opposed to periodic
chains where the local magnetization is expected to exponentially decay to
zero Dağ and Sun (2021). Whether this early time q.s. plateau has a finite
lifetime and the system eventually thermalizes is an interesting and important
question. Because we focus on whether one could utilize the value of the q.s.
regime to probe QCP, an answer to this question is beyond the scope of our
paper. We will characterize the q.s. regime that emerges in open-boundary
nonintegrable TFIC in Secs. IV and V, respectively for short-range and long-
range interacting models.

##  III The integrable TFIC

In this section, we present our analytical and numerical techniques on the
integrable TFIC with $\Delta=0$, which includes (A) the cluster theorem
technique, (B) a proof on the origin of the q.s. regime; (C) the derivation of
a series expression for the edge magnetization; and finally (D) the
numerically extracted single-site nonequilibrium phase diagrams based on the
cluster theorem technique.

###  III.1 Quench dynamics of single-site magnetization via the cluster
theorem

Figure 2: Single-site magnetization $|C_{r}(t)|$ at $r=6$ for different system
sizes (see legend) at $h=0.95$. The magnetization is plotted until
corresponding $t_{l}$, the breakdown time of the cluster theorem, and in this
region the data for different system sizes collapse on each other. This
suggests that we effectively simulate the thermodynamic limit with a finite-
size system at the expense of a breakdown time.

We start by mapping the integrable TFIC to a noninteracting fermionic model in
$1$D via the Jordan–Wigner transformation Sachdev (2001),

$\displaystyle\sigma^{z}_{r}$ $\displaystyle=$
$\displaystyle-\prod_{s<r}\left(1-2c_{s}^{\dagger}c_{s}\right)\left(c_{r}+c_{r}^{\dagger}\right),$
(2) $\displaystyle\sigma^{x}_{r}$ $\displaystyle=$ $\displaystyle
1-2c_{r}^{\dagger}c_{r},$ $\displaystyle\sigma^{y}_{r}$ $\displaystyle=$
$\displaystyle-i\prod_{s<r}\left(1-2c_{s}^{\dagger}c_{s}\right)\left(c_{r}-c_{r}^{\dagger}\right),$

to obtain the noninteracting Hamiltonian

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{r}\big{[}-J\big{(}c^{\dagger}_{r}c_{r+1}+c_{r}^{\dagger}c_{r+1}^{\dagger}+\text{H.c.}\big{)}+2hc_{r}^{\dagger}c_{r}\big{)}\big{]}.$
(3)

Expressing the local magnetization of a bulk spin $\Braket{\sigma_{r}^{z}(t)}$
in this noninteracting picture brings a string of operators, and hence is not
tractable. Instead, one could utilize the Wick theorem and Pfaffian formalism
for two-point correlators to encode the string of operators into an anti-
symmetric matrix Sachdev (2001); Calabrese  _et al._ (2012). This direction of
solution requires rewriting a single-site observable in terms of two-point
correlators Dağ and Sun (2021) and we therefore invoke the cluster theorem
Calabrese  _et al._ (2012),

$\displaystyle\Braket{\sigma_{r}^{z}(t)\sigma_{N-r+1}^{z}(t)}\approx\Braket{\sigma_{r}^{z}(t)}\Braket{\sigma_{N-r+1}^{z}(t)},$
(4)

where $r\ll N/2$ is close to the left boundary. The cluster theorem breaks
down outside of the lightcone, and the corresponding breakdown time $t_{l}$
can be understood in the context of operator spread between two well separated
sites $r$ and $N-r+1$ that eventually get correlated with each other. This
time can be estimated based on the maximum quasiparticle velocities $v_{q}$,
as $t_{l}=\Delta x/(2v_{q})$ where $\Delta x=N-2r+1$ is the distance between
two spins that are equidistant from the symmetry center of an open-boundary
chain in Eq. (4), which is the middle of the chain. For such symmetrically
placed sites, the non-equilibrium response is the same, and one can thus write

$\displaystyle\Braket{\sigma_{r}^{z}(t)}=\sqrt{\Braket{\sigma_{r}^{z}(t)\sigma_{N-r+1}^{z}(t)}}\equiv
C_{r}(t).$ (5)

Therefore, we can extract the dynamical evolution of a spin at site $r$ from
the equal-time two-point correlators of sites $r$ and $N-r+1$. The latter can
be written in terms of auxiliary operators $\phi^{\pm}_{r}=c_{r}^{\dagger}\pm
c_{r}$ as

$\displaystyle\Braket{\sigma_{r}^{z}(t)\sigma_{N-r+1}^{z}(t)}$
$\displaystyle=$
$\displaystyle\bigg{\langle}\phi_{r}^{-}(t)\bigg{(}\prod_{s=r+1}^{N-r}\phi_{s}^{+}(t)\phi_{s}^{-}(t)\bigg{)}$
(6) $\displaystyle\times$ $\displaystyle\phi_{N-r+1}^{+}(t)\bigg{\rangle}.$

This is, in fact, the expectation value of the so-called string order
parameter (SOP) Budich and Heyl (2016). It can be calculated by invoking
Wick’s theorem, which allows one to re-express the above expectation value as
a sum over products of elementary contractions, which in turn is the Pfaffian
of an appropriately constructed antisymmetric matrix $T(t)$ i.e.
$(C_{r}(t))^{2}=\mathrm{Pf}(T(t))$ Sachdev (2001); Calabrese  _et al._ (2012).
Although this is in general a complex number, it follows from Eq. (4) that for
$t<t_{l}$ and the ordered phase we actually have
$(C_{r}(t))^{2}\in\mathbb{R+}$, so that we may compute
$C_{r}(t)=\lvert\mathrm{Pf}(T(t))\rvert^{1/2}=\lvert\sqrt{\mathrm{det}(T(t))}\rvert^{1/2}$.
This is advantageous since it is numerically more efficient to calculate
determinants as compared to Pfaffians. However note that in the disordered
phase and close to QCP in the ordered phase, $C_{r}(t)\in\mathbb{R}$. In these
cases, we compute $|C_{r}(t)|$ with the cluster theorem method. To construct
the matrix $T(t)$ we only need to calculate all possible elementary two-point
contractions $\Braket{\phi_{a}^{p}(t)\phi_{b}^{q}(t)}$ where $p,q=\pm$ and
$r\leq a\leq b\leq N-r+1$. Now we briefly review how to implement a sudden
quench in this picture by following Ref. Uhrich  _et al._ (2020).

Quenching from an initial Hamiltonian $H_{i}$ with transverse field $h_{i}$,
we first solve

$\displaystyle H_{i}=\sum_{k}E_{k}^{i}\alpha_{k}^{\dagger}\alpha_{k},$ (7)

where $E^{i}_{k}$ are the single particle energies, $\alpha_{k}$ and
$\alpha_{k}^{\dagger}$ are the new annihilation and creation operators,
respectively, after diagonalization where the ground state
$\ket{\psi_{0}}_{\alpha}$ satisfies $\alpha_{k}\ket{\psi_{0}}_{\alpha}=0$. The
solution, in a general form, follows as Lieb  _et al._ (1961)

$\displaystyle\left(\begin{array}[]{c}\alpha\\\
\alpha^{\dagger}\end{array}\right)=\left(\begin{array}[]{c c}G_{i}&F_{i}\\\
F_{i}&G_{i}\end{array}\right)\left(\begin{array}[]{c}c_{i}\\\
c^{\dagger}_{i}\end{array}\right),$ (14)

where $c_{i}=(c_{1},c_{2},\cdots,c_{N})^{T}$ and similarly for the creation
operator $c_{i}^{\dagger}$. The expressions for the $G_{i}$ and $F_{i}$ block
matrices will follow shortly. By solving the eigensystem of

$\displaystyle\left[(A_{i}-B_{i})(A_{i}+B_{i})\right]\Ket{\Phi_{k}^{i}}=(E_{k}^{i})^{2}\Ket{\Phi_{k}^{i}},$
(15)

we obtain the eigenenergies $E_{k}^{i}$ and eigenvectors $\Ket{\Phi_{k}^{i}}$.
Note that $A_{i}$ and $B_{i}$ are the nearest neighbor hopping and pairing
terms of the Hamiltonian, respectively. Then the Hamiltonian could be written
as,

$\displaystyle H_{i}=\left(\begin{array}[]{c c}A_{i}&B_{i}\\\
B_{i}^{\dagger}&-A_{i}\end{array}\right),$ (18)

in the $(c\hskip 2.84526ptc^{\dagger})^{T}$ basis. Next we use the eigensystem
$(E_{k}^{i},\Ket{\Phi_{k}^{i}})$ to obtain

$\displaystyle\Ket{\Psi_{k}^{i}}=\frac{1}{E_{k}^{i}}\left[\Bra{\Phi_{k}^{i}}(A_{i}-B_{i})\right]^{T}.$
(19)

Let us emphasize that the expression above would not work for zero edge modes
where $E_{k}^{i}=0$, and handling the case with zero mode is crucial for us.
This is because we are primarily interested in the quenches from and to the
ordered phase in an open-boundary TFIC. In this case, we set
$\Ket{\Psi_{k}^{i}}=-\Ket{\Phi_{k}^{i}}$ Lieb  _et al._ (1961) and make sure
that the resulting state is linearly independent from the other modes. Finally
we calculate the matrices $G_{i}$ and $F_{i}$ in terms of $\Ket{\Phi_{k}^{i}}$
and $\Ket{\Psi_{k}^{i}}$. Defining matrices

$\displaystyle\Phi_{i}$ $\displaystyle=$
$\displaystyle\left[\Ket{\Phi_{1}^{i}}\hskip
2.84526pt\Ket{\Phi_{2}^{i}}\cdots\Ket{\Phi_{N}^{i}}\right],$
$\displaystyle\Psi_{i}$ $\displaystyle=$
$\displaystyle\left[\Ket{\Psi_{1}^{i}}\hskip
2.84526pt\Ket{\Psi_{2}^{i}}\cdots\Ket{\Psi_{N}^{i}}\right],$

the block matrices follow

$\displaystyle G_{i}=\frac{1}{2}\left(\Phi_{i}^{T}+\Psi_{i}^{T}\right),\hskip
2.84526ptF_{i}=\frac{1}{2}\left(\Phi_{i}^{T}-\Psi_{i}^{T}\right).$ (20)

Similarly for the final Hamiltonian $H_{f}$ with transverse field $h$, one
could write

$\displaystyle\left(\begin{array}[]{c}\beta\\\
\beta^{\dagger}\end{array}\right)=\left(\begin{array}[]{c c}G_{f}&F_{f}\\\
F_{f}&G_{f}\end{array}\right)\left(\begin{array}[]{c}c_{f}\\\
c^{\dagger}_{f}\end{array}\right),$ (27)

together with corresponding $\Phi_{f}$ and $\Psi_{f}$. We calculate the
transfer matrices with the help of these block matrices,

$\displaystyle T_{1}=G_{f}G_{i}^{T}+F_{f}F_{i}^{T},$ $\displaystyle
T_{2}=G_{f}F_{i}^{T}+F_{f}G_{i}^{T}.$

These transfer matrices are used to finally calculate the matrix elements
${}_{\alpha}\Bra{\psi_{0}}\left[\phi^{p}_{a}\phi^{q}_{b}\right]_{\beta}\Ket{\psi_{0}}_{\alpha}$
of the antisymmetric matrix $T(t)$, where Greek subscripts denote in which
basis we have the states and the operators. Because we would like to make use
of $\alpha_{k}\Ket{\psi_{0}}_{\alpha}=0$, we use the above transfer matrices
to write $\left[\phi^{p}_{a}\phi^{q}_{b}\right]_{\beta}$ in the $\alpha$ basis
as

$\displaystyle\left[\phi^{\pm}_{b}\right]_{\beta}\Ket{\psi_{0}}_{\alpha}$
$\displaystyle=$ $\displaystyle\left[c_{b}^{\dagger}(t)\pm
c_{b}(t)\right]_{\beta}\Ket{\psi_{0}}_{\alpha},$ $\displaystyle=$
$\displaystyle\left[(G_{f}^{T}\pm F_{f}^{T})\left(e^{i\mathcal{E}t}T_{1}\pm
e^{-i\mathcal{E}t}T_{2}\right)\alpha^{\dagger}\right]_{b}\Ket{\psi_{0}}_{\alpha},$

where $\mathcal{E}$ is a diagonal matrix with eigenenergies of the final
Hamiltonian as the entries, $\mathcal{E}=\text{diag}[E_{1}^{f}\hskip
2.84526ptE_{2}^{f}\cdots E_{N}^{f}]$. Based on this formulation, we construct
matrices $\mathcal{M}_{q}(t)$ explicitly,

$\displaystyle\mathcal{M}_{+}(t)$ $\displaystyle=$
$\displaystyle\Phi_{f}\left(e^{-i\mathcal{E}t}T_{1}+e^{i\mathcal{E}t}T_{2}\right),$
$\displaystyle\mathcal{M}_{-}(t)$ $\displaystyle=$
$\displaystyle\left(T_{1}^{T}e^{i\mathcal{E}t}-T_{2}^{T}e^{-i\mathcal{E}t}\right)\Psi_{f}^{T},$
(28)

to utilize in the following contractions,

$\displaystyle\Braket{\phi^{+}_{a}(t)\phi^{+}_{b}(t)}$ $\displaystyle=$
$\displaystyle[\mathcal{M}_{+}(t)\mathcal{M}_{+}(t)^{\dagger}]_{ab},$
$\displaystyle\Braket{\phi^{-}_{a}(t)\phi^{-}_{b}(t)}$ $\displaystyle=$
$\displaystyle-[\mathcal{M}_{-}^{\dagger}(t)\mathcal{M}_{-}(t)]_{ab},$
$\displaystyle\Braket{\phi^{+}_{a}(t)\phi^{-}_{b}(t)}$ $\displaystyle=$
$\displaystyle[\mathcal{M}_{+}(t)\mathcal{M}_{-}(t)]_{ab},$
$\displaystyle\Braket{\phi^{-}_{a}(t)\phi^{+}_{b}(t)}$ $\displaystyle=$
$\displaystyle-[\mathcal{M}_{-}^{\dagger}(t)\mathcal{M}_{+}^{\dagger}(t)]_{ab}.$
(29)

Now we can construct the antisymmetric matrix $T(t)$ at time $t$ with the
matrix elements $T_{ks}(t)=\Braket{\phi_{a}^{p}(t)\phi_{b}^{q}(t)}$, where
$1\leq k<s\leq 2\Delta x$, $p=+(-)$ for $k$ even(odd) and $q=+(-)$ for $s$
even(odd). The relation between parameters $a,b$ and $k,s$ reads $a=r+\lfloor
k/2\rfloor$ and $b=r+\lfloor s/2\rfloor$, because $r\leq a\leq b\leq N-r+1$.
Having constructed $T(t)$, one can then extract
$C_{r}(t)=|\sqrt{\det{T(t)}}|^{1/2}$, as discussed below Eq. (6). Let us note
that this quench formalism for noninteracting fermions was recently employed
in Refs. Uhrich  _et al._ (2020); Dağ and Sun (2021). We use this method to
obtain the numerical results presented for the critically prethermal regime in
Ref. Dağ  _et al._ (2021).

The breakdown time of the method, $t_{l}$, is also the time when finite-size
effects kick in. Therefore, the method is naturally immune to these finite-
size effects, and the local magnetization for different sizes collapses on
each other, Fig. 2. This is similar in spirit to the separation timescale of
single-site observables which is the time when the quasi-particles reach the
end of the chain Dağ and Sun (2021) and the finite-size effects in this time
interval are exponentially suppressed in system size Wang  _et al._ (2021).

In the next section, we analytically show that (i) the system needs to support
zero modes to stabilize a q.s. regime and (ii) how bulk probe sites are
susceptible to these zero modes which then gives rise to the detection of QCP
at an arbitrary site.

###  III.2 Requirement of zero modes to observe a quasi-stationary regime

In order to show the requirement of zero modes to stabilize a q.s. regime, we
have to adopt an approach different than the cluster theorem for studying the
quench dynamics. This alternative approach, proposed by Ref. Yang (1952),
defines the single-site magnetization as

$\displaystyle C_{r}(t)$ $\displaystyle=$
$\displaystyle(-1)^{r-1}\mathbin{{}_{\alpha}{\Bra{\psi_{0}}}}\left(\prod_{k}^{r-1}\phi_{k}^{+}(t)\phi_{k}^{-}(t)\right)\phi_{r}^{+}(t)\Ket{\psi_{1}}_{\alpha},$

where $\Ket{\psi_{1}}_{\alpha}=\alpha^{\dagger}_{N}\Ket{\psi_{0}}_{\alpha}$ is
the single-particle excited state. The auxiliary operators that are already
defined in Sec. III.1 can be written in the basis of the initial Hamiltonian
in terms of $\alpha_{k}$,

$\displaystyle\phi_{l}^{+}$ $\displaystyle=$
$\displaystyle\sum_{k}^{N}\Phi_{k}^{i}(l)(\alpha_{k}^{\dagger}+\alpha_{k})\equiv\phi_{2l-1},$
(30) $\displaystyle\phi_{l}^{-}$ $\displaystyle=$
$\displaystyle\sum_{k}^{N}\Psi_{k}^{i}(l)(\alpha_{k}^{\dagger}-\alpha_{k})\equiv\phi_{2l}.$
(31)

The parenthesis in $\Phi_{k}^{i}(l)$ and $\Psi_{k}^{i}(l)$ denote the
$l^{\text{th}}$ element of the corresponding vector. The new notation defined
in Eqs. (30) and (31) make the following equations easier to follow,

$\displaystyle C_{1}(t)$ $\displaystyle=$
$\displaystyle\sum_{n,k}\bigg{[}P_{1,2n-1}(t)\Phi_{k}^{i}(n)\mathbin{{}_{\alpha}{\Bra{\psi_{0}}}}(\alpha_{k}^{\dagger}+\alpha_{k})\alpha_{N}^{\dagger}\Ket{\psi_{0}}_{\alpha}$
(32) $\displaystyle+$ $\displaystyle
P_{1,2n}(t)\Psi_{k}^{i}(n)\mathbin{{}_{\alpha}{\Bra{\psi_{0}}}}(\alpha_{k}^{\dagger}-\alpha_{k})\alpha_{N}^{\dagger}\Ket{\psi_{0}}_{\alpha}\bigg{]},$
$\displaystyle=$
$\displaystyle\sum_{n}\left[P_{1,2n-1}(t)\Phi_{N}^{i}(n)-P_{1,2n}(t)\Psi_{N}^{i}(n)\right],$

where $P_{1,2n-1}$ and $P_{1,2n}$ are the single-particle propagators of the
quench Hamiltonian Iglói  _et al._ (2013),

$\displaystyle P_{2l-1,2k-1}(t)$ $\displaystyle=$
$\displaystyle\sum_{q}\cos(E_{q}t)\Phi_{q}^{f}(l)\Phi^{f}_{q}(k),$
$\displaystyle P_{2l-1,2k}(t)$ $\displaystyle=$
$\displaystyle-\sum_{q}\sin(E_{q}t)\Phi^{f}_{q}(l)\Psi^{f}_{q}(k),$
$\displaystyle P_{2l,2k-1}(t)$ $\displaystyle=$
$\displaystyle\sum_{q}\sin(E_{q}t)\Phi^{f}_{q}(k)\Psi^{f}_{q}(l),$
$\displaystyle P_{2l,2k}(t)$ $\displaystyle=$
$\displaystyle\sum_{q}\cos(E_{q}t)\Psi^{f}_{q}(l)\Psi^{f}_{q}(k).$

To be concrete, let us choose the initial state as the ground state of the
Hamiltonian with $h_{i}=0$. This results in

$\displaystyle\Phi_{k}^{i}(n)$ $\displaystyle=$
$\displaystyle\begin{cases}1,\hskip 14.22636pt&\text{if }n=k+1\\\ 0,\hskip
14.22636pt&\text{if }n\neq k+1,\\\ \end{cases}$ (33)
$\displaystyle\Psi_{k}^{i}(n)$ $\displaystyle=$
$\displaystyle\begin{cases}-1,\hskip 14.22636pt&\text{if }n=k\\\ 0,\hskip
14.22636pt&\text{if }n\neq k.\\\ \end{cases}$ (34)

When Eqs. (33) and (34) are substituted into Eq. (32), one obtains

$\displaystyle C_{1}(t)$ $\displaystyle=$ $\displaystyle
P_{1,1}(t)+P_{1,2N}(t),$ $\displaystyle=$
$\displaystyle\sum_{q}\bigg{(}\cos(E_{q}t)\Phi_{q}^{f}(1)\Phi^{f}_{q}(1)$
$\displaystyle-$
$\displaystyle\sin(E_{q}t)\Phi^{f}_{q}(1)\Psi^{f}_{q}(N)\bigg{)}.$ (36)

In the infinite time limit, $t\to\infty$, $E_{q}=0$ is the only nonzero
contribution. Defining the q.s. value $C_{1}(t\rightarrow\infty)\equiv
C_{1}^{qs}(h)$, we obtain

$\displaystyle
C_{1}^{qs}(h)=\sum_{E_{q}=0}\Phi_{q}^{f}(1)\Phi_{q}^{f}(1)=1-h^{2}.$ (37)

This means that the q.s. regime in the edge magnetization is solely due to the
presence of edge modes in the Kitaev chain, Eq. (3). We plot Eq. (37) for a
system size $N=500$ in Fig. 3a, where we show that it coincides perfectly with
$1-h^{2}$ (solid-yellow).

(a)

(b)

Figure 3: (a) Eq. (37) plotted with respect to $h$ (blue-squares) which
coincides with the function $1-h^{2}$ (solid-yellow). (b) Absolute values of
the first (red-solid) and the third term contributions (yellow-dotted) in Eq.
(48) for $N=1500$ spins at $h=0.95$. Their summation (without the absolute
values) gives the entire time evolution for the second spin whose absolute
value is plotted with blue dashed-dotted.

Next, we aim to understand the role of zero modes in the q.s. regime of a bulk
spin, e.g., the second site $r=2$.

$\displaystyle C_{2}(t)$ $\displaystyle=$
$\displaystyle-\mathbin{{}_{\alpha}{\Bra{\psi_{0}}}}\phi_{1}^{+}(t)\phi_{1}^{-}(t)\phi_{2}^{+}(t)\Ket{\psi_{1}}_{\alpha}$
(38) $\displaystyle=$ $\displaystyle\left|\begin{array}[]{c c c
c}0&-\mathcal{P}^{(1,2)}(t)&-\mathcal{P}^{(1,3)}(t)&-\mathcal{P}^{(1)}(t)\\\
\mathcal{P}^{(1,2)}(t)&0&-\mathcal{P}^{(2,3)}(t)&-\mathcal{P}^{(2)}(t)\\\
\mathcal{P}^{(1,3)}(t)&\mathcal{P}^{(2,3)}(t)&0&-\mathcal{P}^{(3)}(t)\\\
\mathcal{P}^{(1)}(t)&\mathcal{P}^{(2)}(t)&\mathcal{P}^{(3)}(t)&0\end{array}\right|^{1/2},$
(43)

where we have used Wick’s theorem, similar to Sec. III.1. The definitions of
this anti-symmetric matrix’s elements read

$\displaystyle\mathcal{P}^{(2r-1)}(t)$ $\displaystyle\equiv$
$\displaystyle\mathbin{{}_{\alpha}{\Bra{\psi_{0}}}}\phi_{r}^{+}(t)\alpha_{N}^{\dagger}\Ket{\psi_{0}}_{\alpha},$
(44) $\displaystyle\mathcal{P}^{(2r)}(t)$ $\displaystyle\equiv$
$\displaystyle\mathbin{{}_{\alpha}{\Bra{\psi_{0}}}}\phi_{r}^{-}(t)\alpha_{N}^{\dagger}\Ket{\psi_{0}}_{\alpha},$
(45) $\displaystyle\mathcal{P}^{(2r-1,2r^{\prime}-1)}(t)$
$\displaystyle\equiv$
$\displaystyle\mathbin{{}_{\alpha}{\Bra{\psi_{0}}}}\phi_{r}^{+}(t)\phi_{r^{\prime}}^{+}(t)\Ket{\psi_{0}}_{\alpha},$
(46) $\displaystyle\mathcal{P}^{(2r-1,2r^{\prime})}(t)$ $\displaystyle\equiv$
$\displaystyle\mathbin{{}_{\alpha}{\Bra{\psi_{0}}}}\phi_{r}^{+}(t)\phi_{r^{\prime}}^{-}(t)\Ket{\psi_{0}}_{\alpha},$
$\displaystyle\mathcal{P}^{(2r,2r^{\prime}-1)}(t)$ $\displaystyle\equiv$
$\displaystyle\mathbin{{}_{\alpha}{\Bra{\psi_{0}}}}\phi_{r}^{-}(t)\phi_{r^{\prime}}^{+}(t)\Ket{\psi_{0}}_{\alpha},$
$\displaystyle\mathcal{P}^{(2r,2r^{\prime})}(t)$ $\displaystyle\equiv$
$\displaystyle\mathbin{{}_{\alpha}{\Bra{\psi_{0}}}}\phi_{r}^{-}(t)\phi_{r^{\prime}}^{-}(t)\Ket{\psi_{0}}_{\alpha}.$

The two-point correlators, Eq. (46) can be calculated with Eq. (29). We know
the form of Eq. (44) from the expression for edge magnetization, Eq. (III.2).
Hence, Eqs. (44) and (45) follow similarly,

$\displaystyle\mathcal{P}^{(2r-1)}(t)$ $\displaystyle=$ $\displaystyle
P_{2r-1,1}(t)+P_{2r-1,2N}(t),$ $\displaystyle\mathcal{P}^{(2r)}(t)$
$\displaystyle=$ $\displaystyle P_{2r,1}(t)+P_{2r,2N}(t).$ (47)

Then, it is straightforward to show

$\displaystyle C_{2}(t)$ $\displaystyle=$
$\displaystyle\mathcal{P}^{(1)}(t)\mathcal{P}^{(2,3)}(t)-\mathcal{P}^{(2)}(t)\mathcal{P}^{(1,3)}(t)$
(48) $\displaystyle+$
$\displaystyle\mathcal{P}^{(3)}(t)\mathcal{P}^{(1,2)}(t).$

In the infinite time limit, $C_{2}(t\rightarrow\infty)$ brings energy
conditions of ${E_{\alpha}\pm E_{\beta}\pm E_{\gamma}=0}$ where all energies
are nonnegative $E_{\alpha}\geq 0$. Then, so long as (i) the single particle
eigenstates are not all product states and (ii) the single particle spectrum
is nondegenerate, all of which is satisfied when $h\neq 0$, the energy
condition holds only when $E_{\alpha}=E_{\beta}=E_{\gamma}=0$. Therefore, to
observe nonzero magnetization in the long time limit, even at bulk probe
sites, the spectrum has to support zero modes. An important difference in the
infinite time limit at $r=2$ from the edge magnetization in Eq. (37) is that
there is a nonzero contribution coming from the second site in the fermionic
chain in addition to the edge. While the contribution of the middle term in
Eq. (48) is always zero because $\Psi_{q|E_{q}=0}^{f}(1)=0$ in
$\mathcal{P}^{(2)}(t)$, the contribution of the first term decreases as the
transverse field increases as opposed to increasing contribution from the
third term. Close to QCP, these two contributions become equal, which is
depicted in Fig. 3b. These observations suggest that the zero mode spreads
across the chain as we approach the QCP, as expected. Importantly, one can
observe the prethermal and q.s. regimes in both terms, shedding light on (i)
why an arbitrary bulk spin is still a good probe site to detect QCP; and (ii)
why a universal collapse is possible in the prethermal regime of different
probe sites Dağ  _et al._ (2021).

Let us note that one could straightforwardly generalize this alternative
numerical method to calculate the local magnetization at all sites. Fig. 1a
shows the application of this method to calculate the local magnetization at
site $r=3$ for $N=1500$ spins. To mathematically show the origin of zero modes
for probe sites $r>2$, one has to add one more condition to the condition list
above, which is the requirement of nondegenerate energy gaps Gogolin and
Eisert (2016).

(a)

(b)

Figure 4: Exact diagonalization results for the coherence time of $C_{3}(t)$
at a system size of $N=14$. Subfigure (a) depicts certain nonintegrable models
at $h=0.5$ (see legend), whereas (b) gives a two-dimensional color plot of the
long-time value of the q.s. regime with respect to external field $h$ and the
interaction strength $\Delta$.

Finally, we discuss the coherence times of the local magnetization at a bulk
site $C_{3}(t)$ (Fig. 4a), with respect to time for different interaction
strengths $\Delta$. Then we map the q.s. values at a field strength $h$ and
the interaction strength $\Delta$ in Fig. 4b. The behavior is monotonous
everywhere between $0<h<1$ and $0<\Delta<2$. An important evidence of strong
zero modes Kemp  _et al._ (2017) is the presence of resonances due to the
competition between n.n. and n.n.n. terms in Eq. (1), which would result in a
non-monotonous trend of the steady state value with respect to $\Delta$. The
absence of such behavior seen in Fig. 4, despite our analytical proof earlier
on the necessity of zero modes to stabilize a q.s. regime, is likely a result
of the quench protocol. This is because, in a typical setup where the effects
of strong zero modes are demonstrated, e.g., infinitely long or finite but
long coherence times, the entire many-body spectrum is relevant due to
infinite-temperature initial state Kemp  _et al._ (2017); Gómez-Ruiz  _et al._
(2018); Dağ  _et al._ (2020). Although a quench from a polarized state to an
arbitrary $h$ in the ordered phase still generates excitations to higher
energy levels, Fig. 4 clearly demonstrates that only a part of the many-body
spectrum is actually relevant in the observation of the q.s. regime.

###  III.3 Derivation of the series expression for the edge magnetization

In order to write the edge magnetization, Eq. (36), in terms of the physical
parameters, e.g., transverse field $h$, one has to analytically express
$\Phi_{q}^{f}$ and $\Psi_{q}^{f}$. In this section, we will follow an
alternative route to derive an analytical expression for $C_{1}(t)$ in the
thermodynamic limit. First we rewrite the Hamiltonian of Eq. (3) at $\Delta=0$
in terms of auxilary operators $\phi_{r}^{\pm}$,

$\displaystyle
H=\sum_{r}(-J\phi_{r}^{-}\phi_{r+1}^{+}+h\phi_{r}^{+}\phi_{r}^{-}).$ (49)

Then we use the Baker–Campbell–Hausdorff formula

$\displaystyle e^{sA}Be^{-sA}=B+s[A,B]+\frac{s^{2}}{2!}[A,[A,B]]+\dots,$ (50)

to calculate $\phi_{1}^{+}(t)=e^{iHt}\phi_{1}^{+}e^{-iHt}$. The following
commutators are obtained,

$\displaystyle[H,\phi_{1}^{+}]$ $\displaystyle=$ $\displaystyle
2h\phi_{1}^{-},$ $\displaystyle\left[H,\phi_{2}^{+}\right]$ $\displaystyle=$
$\displaystyle 2(J\phi_{1}^{-}+h\phi_{2}^{-}),$
$\displaystyle\left[H,\phi_{3}^{+}\right]$ $\displaystyle=$ $\displaystyle
2(J\phi_{2}^{-}+h\phi_{3}^{-}),$ $\displaystyle\dots$
$\displaystyle\left[H,\phi_{1}^{-}\right]$ $\displaystyle=$
$\displaystyle-2(h\phi_{1}^{+}+J\phi_{2}^{+}),$
$\displaystyle\left[H,\phi_{2}^{-}\right]$ $\displaystyle=$
$\displaystyle-2(h\phi_{2}^{+}+J\phi_{3}^{+}),$ $\displaystyle\dots$

and used to calculate the terms of Eq. (50),

$\displaystyle[H,[H,\phi_{1}^{+}]]$ $\displaystyle=$
$\displaystyle-4(h^{2}\phi_{1}^{+}+Jh\phi_{2}^{+}),$
$\displaystyle\left[H,\left[H,\left[H,\gamma_{1}\right]\right]\right]$
$\displaystyle=$
$\displaystyle-8[(h^{3}+J^{2}h)\phi_{1}^{-}+Jh^{2}\phi_{2}^{-}],$
$\displaystyle\left[H,\left[H,\left[H,\left[H,\gamma_{1}\right]\right]\right]\right]$
$\displaystyle=$ $\displaystyle 16[(h^{4}+J^{2}h^{2})\phi_{1}^{+}+(hJ^{3}$
$\displaystyle+$ $\displaystyle 2h^{3}J)\phi_{2}^{+}+J^{2}h^{2}\phi_{3}^{+}],$
$\displaystyle\dots$ (51)

Since we calculate
$C_{1}(t)=\mathbin{{}_{\alpha}{\Bra{\psi_{0}}}}\phi_{1}^{+}(t)\Ket{\psi_{1}}_{\alpha}$
quenched from a polarized state, the initial state dictates
$\Braket{\phi_{1}^{+}}=1$, $\Braket{\phi_{N}^{-}}=-i$ and the remaining
expectation values are zero. Hence, Eq. (51) read,

$\displaystyle\Braket{[H,[H,\gamma_{1}]]}_{\alpha}$ $\displaystyle=$
$\displaystyle-2^{2}h^{2},$
$\displaystyle\Braket{\left[H,\left[H,\left[H,\gamma_{1}\right]\right]\right]}_{\alpha}$
$\displaystyle=$ $\displaystyle 0$
$\displaystyle\Braket{\left[H,\left[H,\left[H,\left[H,\gamma_{1}\right]\right]\right]\right]}_{\alpha}$
$\displaystyle=$ $\displaystyle 2^{4}(h^{4}+J^{2}h^{2})$
$\displaystyle\Braket{\left[H,\left[H,\left[H,\left[H,\left[H,\gamma_{1}\right]\right]\right]\right]\right]}_{\alpha}$
$\displaystyle=$ $\displaystyle 0$
$\displaystyle\Braket{\left[H,\left[H,\left[H,\left[H,\left[H,\left[H,\gamma_{1}\right]\right]\right]\right]\right]\right]}_{\alpha}$
$\displaystyle=$ $\displaystyle-2^{6}(h^{6}$ $\displaystyle+$ $\displaystyle
3h^{4}J^{2}+h^{2}J^{4})$ $\displaystyle\dots$

resulting in the following series solution

$\displaystyle
C_{1}(t)=\sum_{m=0}\frac{(-1)^{m}}{(2m)!}(2t)^{2m}N_{m}(h^{2}),$ (52)

where $N_{m}(x)$ are called the Narayana polynomials

$\displaystyle N_{m}(h)=\sum_{n=1}^{m}N_{mn}h^{n},\hskip
8.53581ptN_{mn}=\frac{1}{m}\binom{m}{n-1}\binom{m}{n}.$

Let us note that one obtains exactly the same series solution in the
calculation of two-time correlators of the edge spin at infinite-temperature
Gómez-Ruiz  _et al._ (2018).

(a)

(b)

Figure 5: The local magnetizations of the integrable TFIC (a) at QCP $h=1$ and
(b) in the disordered phase $h=1.2$ for different probe sites (see legend).
All decays as a power-law in long time with $t^{-3/2}$.

At the critical point $h=1$, the series solution Eq. (52) can be written in a
closed form,

$\displaystyle C_{1}(t,h=1)=\frac{J_{1}(4t)}{2t},$ (53)

where $J_{1}(4t)$ is the Bessel function of the first kind. The long-time
asymptotics of $J_{1}(4t)/2t$ is $\propto t^{-3/2}$ which is exactly what we
observe in the numerical calculations, Fig. 5a, not only for the edge
magnetization but also for local magnetization at different probe sites.

Since the Narayana polynomials satisfy the property
$N_{m}(x^{-1})=x^{-m-1}N_{m}(x)$, $C_{1}(t,h)$ satisfies

$C_{1}(t,h)=1-h^{2}+h^{2}C_{1}(ht,h^{-1}),$ (54)

indicating that the edge magnetization behaves in one phase (transverse-field
strength $h$) as it would in the complementary phase (transverse-field
strength $h^{-1}$) but with a scaled time $ht$.

In Ref. Dağ  _et al._ (2021), some of us proposed a general form for the
critical response function $\delta C_{r}(t,h_{n})\equiv
C_{r}(t,h_{n})-C_{r}(t,h_{n}=0)$ where $h_{n}\equiv(h_{c}-h)/h_{c}$,

$\displaystyle\delta
C_{r}(t,h_{n})=C^{qs}_{r}(|h_{n}|)f_{\Delta,h_{i}}(h_{n}t)$ (55)

which holds near $h_{n}=0$ when $t\gg 1$, and the function $f_{0,0}(h_{n}t)$
has been derived based on the edge magnetization in Eq. (52),

$f_{0,0}(h_{n}t)=\frac{1}{2}+\frac{(h_{n}t)}{2}{}_{1}F_{2}\left[\left\\{\frac{1}{2}\right\\};\left\\{\frac{3}{2},2\right\\};-(h_{n}t)^{2}\right]$
(56)

where
${}_{1}F_{2}\left(\left\\{a_{1}\right\\};\left\\{b_{1},b_{2}\right\\};z\right)$
is a generalized Hypergeometric function. Even though Ref. Dağ  _et al._
(2021) focuses on dynamics when quenching to the ferromagnetic phase, i.e.,
$h<1\ (h_{n}>0)$, the function $f_{0,0}(h_{n}t)$ given by Eq. (56) is a good
approximation of $\delta C_{r}(t,h_{n})$ in both phases near QCP,
$|h_{n}|\rightarrow 0$. This suggests that there is a corresponding prethermal
temporal regime in the disordered phase too where dynamics critically slow
down.

Figure 6: The rescaled critical response $\tilde{h}^{-1}\delta C_{1}(t,h)$ is
plotted with respect to the rescaled time $\tilde{h}t$. The blue and red solid
lines show the critical response for the ordered and disordered phases,
respectively. Dashed lines are added to guide the eye.  Figure 7:
Nonequilibrium phase diagram for $\Delta=0$ (a) for different sites
$r=1,3,6,9,12$ at $N=1440$ (see legend). Solid-black line is the analytic
result for edge magnetization in the thermodynamic limit. (b) The same plot
with (a) in semi-logarithmic scale, exhibiting a kink behavior at the QCP
which is nonanalytic in the thermodynamic limit. The small but finite values
at the QCP and in the disordered phase are due to inevitable finite-time
cutoff. These values are expected to vanish when $t_{l}\rightarrow\infty$, see
Fig. 5. (c) For different system sizes, $N=500$ – $2000$, at site $r=6$.
$\overline{C_{r}(t)}$ at $h_{c}$ scales with the system size as $1/N$ for all
sites $r=3,6,9,12$ (see Appendix A).

In Fig. 6, we plot the function $f_{0,0}(h_{n}t)$ as a function of $h_{n}t$.
Since the second term on the right hand side of Eq. (56) is an odd function in
$h_{n}t$, $f_{0,0}(h_{n}t)$ is symmetric about the line $\frac{1}{2}$ in Fig.
6 both for $h_{n}>0$ and $h_{n}<0$. We observe $\delta C_{1}(t=0,h)=h_{n}/2$
in both phases when $t\gg 1$ and $h_{n}t\ll 1$ hold, hence the prethermal
dynamics starts at value $h_{n}/2$ in both phases. $\delta
C_{1}(\infty,h<1)=1-h^{2}$, while $\delta C_{1}(\infty,h>1)=0$, indicating
that the edge magnetization reaches a q.s. value $1-h^{2}$ in the ordered
phase and $0$ in the disordered phase as $t\rightarrow\infty$. Therefore, the
critical response shows a non-analyticity at QCP in the long time limit. We
will harness this non-analyticity in the next sections in the detection of QCP
with single-site magnetizations.

###  III.4 Detection of QCP

To probe the q.s. regime, we study the the time-averaged single-site
magnetization,

$\displaystyle\overline{|C_{r}(t,h)|}=\frac{1}{(t_{l}-t^{*})}\int_{t^{*}}^{t_{l}}dt\,|C_{r}(t,h)|\equiv\bar{C}_{r}(h)$
(57)

where $t^{*}$ is the ultraviolet (short-time, short-distance) cutoff and
$t_{l}$ is the (long-time, long-distance) infrared cutoff Dağ and Sun (2021).
For the numerical results presented in this section, $t_{l}$ is the evolution
time at which the cluster theorem, Sec. III.1, breaks down, and hence as
explained in the previous section $t_{l}=\Delta x/(2v_{q})$. This means that
the infrared cutoff $t_{l}$ is parametric in the transverse field $t_{l}(h)$,
because in the integrable TFIC $v_{q}=2h$ for $h\leq h_{c}$ and $v_{q}=2h_{c}$
for $h>h_{c}$ Calabrese  _et al._ (2012).

For quenches sufficiently far away from the vicinity of the QCP,
$\overline{|C_{r}(t)|}$ matches the q.s. value, as then there is no prethermal
regime Dağ  _et al._ (2021). Each single-site observable in the TFIC
equilibrates around a different value in the q.s. regime. This can be seen in
Fig. 7a that depicts $\overline{|C_{r}(t)|}$ for $r=1-12$, all of which have a
different functional form of $h$. The analytic expression for the value of the
q.s. regime at $r=1$, $C^{qs}_{1}(h)=1-h^{2}$ (Sec. III.3) matches the
corresponding numerical result in Fig. 7a. Fig. 7b shows the same plot in
logarithmic scale that exhibits a kink in $\bar{C}_{r}(h_{c})$ regardless of
the value of $r$.

Whether there is a nonanalytic behavior at QCP for all $r$ could be answered
through a finite-size analysis. In Fig. 7c, we observe that the dynamical
order of $\bar{C}_{6}(h)$ is persistent for $h<h_{c}$ as system size increases
from $N=500$ to $2000$, i.e., its value stays the same with increasing system
size (see Appendix A for other $r$). Although the local order profiles for
each $r$ away from the QCP are different, Fig. 7a, they all approach to QCP
linearly in $h_{n}$ in the thermodynamic limit Dağ  _et al._ (2021). For
$h\geq h_{c}$ the dynamic order vanishes for all $r$ in Fig. 7c, i.e., it has
a decreasing trend with increasing system size. In fact, we can prove that the
local magnetization should vanish in the thermodynamic limit, based on the
simple observation of $t^{-3/2}$ decay of local magnetization for all $r$ both
at the QCP, Fig. 5a and in the disordered phase, Fig. 5b.

$\displaystyle\overline{C(t,h\geq
h_{c})}=\frac{8}{\sqrt{t^{*}}}\left(\frac{N}{h_{c}}\right)^{-1}+\mathcal{O}(N^{-3/2}).$
(58)

The derivation of Eq. (58) and the accompanying numerics of the scaling $1/N$
at the QCP can be found in Appendix A. Therefore, the nonanalyticity at
$h_{c}$—the hallmark of a phase transition—is captured by all single-site
observables $r\ll N/2$. Consequently, we have demonstrated a DPT for different
sites $r\ll N/2$ that reflects the underlying QPT.

##  IV The nonintegrable TFIC with next nearest neighbor coupling

In this section, we first lay out the qMFT formalism that has also been
applied to obtain the results of weak perturbations in Ref. Dağ  _et al._
(2021). Subsequently, we present the $t$-DMRG calculations of both weakly and
strongly interacting TFIC.

###  IV.1 Quench Mean-field theory (qMFT) analysis

(a)  Figure 8: qMFT single-site magnetization of the near-integrable model
with $\Delta=0.1$ (a) for $h\leq 1.143$, an interval of $h$ that is
representative of a dynamically-ordered phase and (b) for $1.146\leq h\leq
1.5$, an interval of $h$ that is representative of a dynamically-disordered
phase at a system size of $N=1440$ at $r=6$. The recurrence attempts observed
here imply the breakdown of the cluster theorem. qMFT single-site
magnetization of the near-integrable model for probe sites $r=1,3,6,9,12$ (c)
at DCP $h_{dc}=1.1437$, (d) in the disordered phase sharing the same legend.

In order to incorporate the qMFT formalism to our quench dynamics that is
described in Sec. III.1, we need to express the n.n.n. term in the TFIC
Hamiltonian in the fermionic picture. This expression reads,

n.n.n. $\displaystyle=$
$\displaystyle\Delta\sum_{r}\left(c_{r}^{\dagger}-c_{r}\right)\left(1-2c_{r+1}^{\dagger}c_{r+1}\right)\left(c_{r+2}+c_{r+2}^{\dagger}\right),$
(59) $\displaystyle=$
$\displaystyle\Delta\sum_{r}\phi^{-}_{r}\phi^{+}_{r+1}\phi^{-}_{r+1}\phi^{+}_{r+2},$

where $\Delta>0$ and $\phi_{r}^{\pm}$ stand for the auxiliary fermions.

In Hartree–Fock expansion, we assume $|\Delta|\ll|J|$, and write Eq. (59) as

$\displaystyle=$
$\displaystyle\Delta\sum_{r}\bigg{[}\Braket{\phi^{-}_{r}\phi^{+}_{r+1}}_{t\rightarrow\infty}\phi^{-}_{r+1}\phi^{+}_{r+2}$
(60) $\displaystyle+$
$\displaystyle\phi^{-}_{r}\phi^{+}_{r+1}\Braket{\phi^{-}_{r+1}\phi^{+}_{r+2}}_{t\rightarrow\infty}-\phi^{-}_{r}\Braket{\phi^{+}_{r+1}\phi^{-}_{r+1}}_{t\rightarrow\infty}\phi^{+}_{r+2}$
$\displaystyle-$
$\displaystyle\Braket{\phi^{-}_{r}\phi^{+}_{r+2}}_{t\rightarrow\infty}\phi^{+}_{r+1}\phi^{-}_{r+1}\bigg{]}.$

Here the $\Braket{\cdot}_{t\rightarrow\infty}$ means that we calculate the
free fermion problem and obtain the correlators with respect to the state in
the q.s. regime when $t\rightarrow\infty$, instead of the ground state which
would be used to calculate the equilibrium QCP. An analytical qMFT formalism
was first introduced in Ref. Titum  _et al._ (2019) to calculate a dynamical
order parameter, based on two-point correlators in a periodic chain, by
utilizing the momentum space representation where the infinite time limit can
indeed be taken. In our numerics for an open-boundary chain, we treat the
largest time point allowed by the cluster theorem, $t_{l}=\Delta x/(2v_{q})$,
as the asymptotically infinite time limit. Note that for an open-boundary
chain one needs to carefully take the edges of the chain into account when
calculating Eq. (60). Using the above expansion, we obtain an effective mean
field Hamiltonian which has slightly stronger n.n. coupling compared to the
free problem, as well as new n.n.n couplings. Further, the effective chemical
potential slightly decreases, which is reasonable when we think about how the
critical point shifts to favor order, e.g., for $\Delta=0.1$, $h_{c}=1.1631\pm
0.0037$ (see Appendix C).

When applying the cluster theorem to the qMFT of the weakly-interacting
nonintegrable TFIC, one needs to estimate the lightcone (correlation) velocity
$v_{q}$ of the model. While for the integrable TFIC this velocity is
analytically known, this is not true when we introduce nonintegrability to the
model. Here we approximate a quasiparticle velocity based on the analytical
prediction of the integrable TFIC: $v_{q}=2h$ for $h\leq h_{c}$ and
$v_{q}=2h_{c}$ for $h>h_{c}$ Calabrese  _et al._ (2012). Since this is only an
approximation, we sometimes exceed the time when the cluster theorem really
breaks down. This means that the distant sites of the chain had already become
correlated with one another. This time can be observed with a recurrence
attempt in the results, e.g., Fig. 8a, which is also a sign of finite-size
effects.

Figures 8a and 8b show single-site magnetization of the near-integrable model
with $\Delta=0.1$ at probe site $r=6$ and for different $h$ across the QCP,
$h_{c}$. One could easily notice the qualitative difference between single-
site magnetizations in different phases: While there is a q.s. regime present
in Fig. 8a for $h$ sufficiently away from the QCP signaling a dynamically-
ordered phase, the magnetization decays as a power-law in Fig. 8b for $h\geq
h_{c}$ indicating a dynamically-disordered phase.

###  IV.2 Detection of QCP in the near integrable TFIC through qMFT analysis

Now we study the single-site magnetization of the near integrable TFIC where
$\Delta=0.1$ in Eq. (1) treated with qMFT and when the system is quenched to
the vicinity of QCP. This model has a QCP at $h_{c}=1.1631\pm 0.0037$ which is
calculated with DMRG (Appendix C).

Figure 9: qMFT nonequilibrium phase diagram of the near-integrable TFIC with
$\Delta=0.1$ (a) for different sites $r=1,3,6,9,12$ in descending order at
$N=1440$ exhibiting a kink at the DCP and solid-black line is a fit function
for the edge magnetization; (b) for different system sizes $N=96$ – $1440$ at
site $r=6$.

Local order profiles follow similarly to the case of integrable TFIC
($\Delta=0$) except for a shift in the QCP to favor order, $h_{c}>1$, as
expected (Fig. 9a). We focus on $\bar{C}_{6}(h)$ in Fig. 9b where we find a
kink at $h_{dc}=1.1437\pm 0.0001<h_{c}$ whose magnetization decreases with
increasing system size $\bar{C}_{r}(h_{c})\propto N^{-1}$. $N^{-1}$ scaling
follows from the fact that the single-site magnetization at $h_{dc}$ decays as
$\propto t^{-3/2}$ (see Fig. 8c and Appendix C). The magnetization in the
disordered phase decays in the same way in Fig. 8d. Consistently, this
behavior is translated to the local order profiles as $\bar{C}_{r}(h\geq
h_{dc})\propto N^{-1}$ suggesting $\bar{C}_{r}(h\geq h_{dc})=0$ in the
thermodynamic limit. On the contrary for $h<h_{dc}$, the dynamical order is
persistent, i.e., $\bar{C}_{r}(h<h_{dc})$ value stays the same with increasing
system size.

The q.s. value of the edge magnetization, which our numerical results access
for quenches far away from the vicinity of the transition, can be fitted well
with a functional form that is reminiscent of that of the integrable TFIC,
$C^{qs}_{1}(h)=\alpha(h_{dc}^{\nu}-h^{\nu})$ for $h\leq h_{dc}$ and zero
otherwise. Furthermore, because a critically prethermal regime with self-
similar dynamics is found also for the same model in Dağ  _et al._ (2021), all
local order profiles in the ordered phase approach to QCP linearly in $h_{n}$
in the thermodynamic limit. Therefore, we observe a non-analytic behavior for
the near-integrable TFIC when it is treated with qMFT. Subsequently, we call
the location of this non-analyticity as a DCP.

We also apply energy gap analysis to the effective qMFT Hamiltonian with
n.n.n. terms in Eq. (59) and find the ground state gap of this effective
Hamiltonian closes at $h_{c}^{\text{qmft}}=1.1438\pm 0.0011$ (see Appendix C).
Let us note that $h_{c}^{\text{qmft}}\sim h_{dc}<h_{c}$, suggesting that the
DCP of the near-integrable model traces back to the QCP of the effective qMFT
Hamiltonian. Therefore, we conclude that the small shift in $h_{dc}$ from
$h_{c}$ likely originates from the qMFT method.

It is important to note that although the qMFT treatment gives rise to a q.s.
regime for long intervals of time for quenches away from the vicinity of the
DCP (Fig. 8a), this is not conclusive evidence for infinitely long-lived
nonthermal behavior in a nonintegrable model. Indeed, MFT is not expected to
adequately capture thermalization as it may neglect fluctuations that are
essential for the latter. However as we will detail in the next subsection,
$t-$DMRG confirms the presence of a transition in the thermodynamic limit.
Since the relaxation time of the system observables critically diverges as we
approach the QCP Dağ  _et al._ (2021), it is likely that an infrared cutoff at
$t_{l}$ might not accurately capture q.s. state in the infinite-limit near
QCP. This could explain why our observed DCP at $h_{dc}=1.1437\pm 0.0001$ for
$\Delta=0.1$ is slightly smaller than the DCP found in Ref. Titum  _et al._
(2019), which is $1.15$.

After determining the DCP, we obtain $\alpha=0.81$ and $\nu=1.81$ for the fit
function of the edge magnetization. This fit function is plotted in Fig. 9a as
a black-solid line.

###  IV.3 $t$-DMRG calculations

(a)

(b)

Figure 10: Order profiles calculated with $t$-DMRG for weakly nonintegrable
TFIC with $\Delta=0.1J$, (a) for different system sizes $N=48,72,96$ at site
$r=3$ when a (infrared) temporal cutoff of $t_{l}=N/2$ is applied. The order
profile is tested against different infrared and ultraviolet temporal cutoffs
and shown to be robust (see text). (b) Finite-time scaling analysis at a
system size of $N=96$ for different infrared cutoffs, shown in the legend,
with a fixed ultraviolet cutoff of $t^{*}=0$. The upward and downward arrows
highlight $h$ points where the order grows or diminishes, respectively with
increasing (a) system size and (b) simulation time.

(a)

(b)

(c)

Figure 11: Order profiles calculated with $t$-DMRG for weakly nonintegrable
TFIC with $\Delta=0.1J$, (a) for different system sizes $N=72,96,120$ at site
$r=6$ when a (infrared) temporal cutoff of $t_{l}=N/2$ is applied. The order
profile is tested against different infrared and ultraviolet temporal cutoffs
and shown to be robust (see text). (b) Finite-time scaling analysis at a
system size of $N=120$ for different infrared cutoffs, shown in the legend,
with a fixed ultraviolet cutoff of $t^{*}=0$. The upward and downward arrows
highlight the $h/J$ points where the order grows or diminishes, respectively,
with increasing (a) system size and (b) simulation time. (c) Zoom over the
phase diagram based on the qMFT result with arrows demonstrating how the phase
diagram of the smallest system size $N=96$ crosses through those of the larger
system sizes before the found DCP $h_{dc}\approx 1.1437$ with the qMFT
analysis.

(a)

(b)

Figure 12: Order profiles calculated with $t$-DMRG for strongly nonintegrable
TFIC with $\Delta=1$ for different system sizes $N=42,48$ at sites (a) $r=3$
and (b) $r=6$ when an infrared temporal cutoff of $t_{l}=N/3$ and ultraviolet
temporal cutoff of $t^{*}=0$ are applied.

(a)

(b)

(c)

Figure 13: Order profiles calculated with $t$-DMRG for strongly nonintegrable
TFIC with $\Delta=1$ based on finite-time analysis at system size $N=48$ for
(a) site $r=3$ and (b) site $r=6$. The infrared cutoffs are shown in the
legends and the ultraviolet cutoff is $t^{*}=0$. (c) Same parameters as (a),
however the maximum bond dimension is set to $\chi_{m}=200$, instead of
$\chi_{m}=100$.

We utilize the ITensor environment ITe  to construct matrix product states
(MPS) and Trotter decomposition for the time evolution of the MPS. We set a
maximum bond dimension $\chi_{m}$ for the resulting compressed MPS, and set
the initial truncation error cutoff for the compression of the MPS as
$\epsilon\approx 10^{-8}$. The truncation error cutoff is adaptive: As the
maximum bond dimension is reached for the resulting MPS, the error cutoff
increases systematically up until a hard error threshold of $\epsilon\approx
10^{-5}$ to be able to access longer times. Setting a maximum allowed bond
dimension thus introduces an error which grows with time. Consequently, we are
confined to early times for which the above interval of the error thresholds
is satisfied.

We study two nonintegrable models with $t$-DMRG: (i) Weakly interacting near-
integrable model, $\Delta=0.1$, whose qMFT results are presented in the
previous subsection and compared to the results obtained from the exact
numerics in this subsection. (ii) A strongly nonintegrable model, $\Delta=1$
where qMFT is inadequate, and we therefore employ $t$-DMRG only. The latter is
numerically exact, but within a given fidelity threshold the accessible
evolution times are limited and naturally shorter than those achievable for
the (near-)integrable model.

Figs. 10a and 10b show the finite-size and finite-time scaling analysis,
respectively for the weakly interacting TFIC at $r=3$. In the rest of this
subsection, we use infrared cutoffs either parametric in the system size
$t_{l}(N)$ (in finite-size analysis) or set to a fixed time $t_{l}$ (in
finite-time analysis).

In Fig. 10a, we observe a dynamical order profile and a crossing point for
system sizes between $N=48$ and $N=96$. The crossing point is found to be in
the interval of $h_{dc}/J\in(1.05,1.1)$ where the dynamical order builds up
for $h/J\leq 1.05$ as the system size increases and it vanishes for $h/J\geq
1.1$. We use $t^{*}=10$ and $t_{l}=N/2$ as ultraviolet and infrared cutoffs in
these figures, however the phase diagram is robust to changes in the temporal
cutoffs, as we tested with values in the intervals $t^{*}\in[8,12]$ and
$t_{l}\in[N/3,N/2]$. In Fig. 10b, we employ a form of finite-time scaling
analysis Titum  _et al._ (2019) where we see a crossing point between order
profiles with different infrared cutoffs at a fixed ultraviolet cutoff
$t^{*}=0$. The crossing point resides in the interval of
$h_{dc}/J\in(1.05,1.1)$ agreeing with the finite-size analysis. For $h/J\leq
1.05$ and $h/J\geq 1.1$, the dynamical order grows or diminishes with longer
simulation times, respectively. Given that the system sizes are constrained in
$t$-DMRG compared to qMFT analysis, we observe a smaller DCP than what we have
found with qMFT analysis. We repeat the same calculation for $r=6$, however
with slightly larger system sizes. Figs. 11a and 11b demonstrate the finite-
size and finite-time scaling analysis, respectively for this parameter set. In
both cases, we find the crossing interval to be $h_{dc}/J\in(1.1,1.15)$.
Finally, we zoom on the nonequilibrium phase diagram of the qMFT analysis in
Fig. 11c to demonstrate that the crossing is present in this method too,
however for small system sizes, e.g., $N=96$ shown with black arrows around
$h/J\approx 1.14$. On the contrary, all system sizes $N\geq 480$ collapse on
the qMFT-predicted DCP suggesting that these finite but large system sizes
effectively simulate the thermodynamic limit. Therefore, we conclude that the
$t$-DMRG method of the simulated system sizes $N\leq 120$ supports the
presence of a DCP, verifying the results of the qMFT method presented
previously.

Next, we consider the strongly nonintegrable TFIC with $\Delta=1$, whose QCP
lies at $h_{c}\approx 2.46$ Dağ and Sun (2021). Figure 12a and 12b show the
local order profiles for $r=3$ and $r=6$, respectively. As the system becomes
strongly interacting, the maximum system size and the evolution times become
more limited. In this set of calculations, we study $N=42$ and $N=48$ with
maximum evolution times of $t\approx 20$. Hence, we apply an infrared cutoff
of $t_{l}=N/3$ and set an ultraviolet (UV) cutoff $t^{*}=0$ to increase the
temporal range of data to average over. Both figures exhibit a crossing point
similar to that of the discussion in the previous paragraph, in the interval
of $h_{dc}/J\in(2.3,2.35)<h_{c}$ where for $h/J<2.3$ the dynamical order
increase with increasing system size, and for $h/J>2.35$ decreases. Therefore,
we conclude that (i) the crossing point seems to be independent of the
measurement site, and (ii) the presence of a crossing point in $t$-DMRG data
hints at the presence of a DCP. We note that more data with larger system
sizes and longer simulations times are required to test the robustness of
these results, determine the exact location of the DCP, and further look for a
critically prethermal regime and dynamical scaling in this strongly
nonintegrable model. Nevertheless, single-site observables close to an edge
seem to be a probe of criticality generally in short-range models, not limited
to noninteracting or weakly-interacting systems.

Finally, we apply the finite-time scaling analysis on the strongly
nonintegrable TFIC with $\Delta=1$ shown in Figs. 13a and 13b for $r=3$ and
$r=6$, respectively. Although there is a well-defined finite-time crossing
point for this set of parameters, the crossing suggests $h_{dc}/J\approx 2.2$
for $r=3$ and $h_{dc}/J\in(2.25,2.3)$ for $r=6$ all of which is less than the
result determined by the finite-size analysis. To check whether the results
depend on the maximum bond dimension, we set $\chi_{m}=200$ for $r=3$, repeat
the calculation on the relevant region in Fig. 13c, and observe that the
crossing point does not change. Although more data is necessary for conclusive
results, we note that the mismatch between finite-size and finite-time
analyses as well as the discrepancy between different sites in finite-time
analysis might point to a change in the light cone structure, e.g., from
linear to power-law, as the transverse field increases and approaches the
equilibrium QPT in finite-size systems.

##  V The quasi-stationary regime in long-range interacting nonintegrable TFIC

The q.s. temporal regime also emerges in the long-range hard-boundary TFIC
with power-law decaying interactions, as was previously noted in Ref.
Neyenhuis  _et al._ (2017) in the context of prethermalization. The
Hamiltonian for the one-dimensional long-range TFIC reads,

$\displaystyle
H=-\sum_{r,r^{\prime}}J(r,r^{\prime})\sigma_{r}^{z}\sigma_{r^{\prime}}^{z}+h\sum_{r}\sigma_{r}^{x},$
(61)

where $J(r,r^{\prime})=J/|r-r^{\prime}|^{\alpha}$. In the limit where
$\alpha=0$, the model becomes integrable with all-to-all interactions, e.g.
LMG model; whereas in the limit of $\alpha\rightarrow\infty$ the model reduces
to short-range n.n. TFIC. For $\alpha\geq 3$, the model belongs to the short-
range Ising universality class Dutta and Bhattacharjee (2001); Knap  _et al._
(2013), and it has algebraic light cones for $\alpha>2$ that approach linear
cones as $\alpha\rightarrow\infty$ Foss-Feig  _et al._ (2015). In this
section, we provide numerical evidence on the presence of a q.s. regime in
TFIC with various $\alpha$, and determine when the q.s. regime breaks down.
Furthermore, we calculate the local order profiles at $\alpha=2.5$, which is a
long-range model with algebraic light cones.

All data presented in this section is obtained from the time-dependent DMRG
method Schollwöck (2005, 2011); Daley  _et al._ (2004); Gobert  _et al._
(2005); mpt  with Krylov time evolution kry . We find convergence for a time-
step of $\Delta t=0.01J$ and a fidelity threshold between $10^{-8}$ –
$10^{-6}$ for our most stringent calculations. Power-law profiles are
approximated as a sum of five exponentials fitted over the length of the
chain, as is optimal for a matrix product operator formulation Crosswhite  _et
al._ (2008).

(a)

(b)

(c)

(d)

Figure 14: Single-site magnetization for $r=N/2$ (spin in the middle of the
chain), $r=6$ (spin close to the boundary) and the total magnetization, when
(a) $\alpha=4$, (b) $\alpha=3$, (c) $\alpha=2.5$ and (d) $\alpha=2$ at
transverse field $h=0.5\tilde{h}_{c}$.

Fig. 14a shows the edge magnetization of single-site observables near the
edge, and in the middle of the chain, as well as the total magnetization
evolving under Hamiltonian Eq. (61) with $\alpha=4$ and $h=0.5\tilde{h}_{c}$
quenched from a polarized state. The unit notation $\tilde{h}_{c}$ will be
explained shortly. We observe that the single-site observable close to the
boundary (red) develops a q.s. regime whereas the observable in the middle of
the chain (blue) exhibits decay. This behavior, consistently, is similar to
what we have observed in the n.n. TFIC model, c.f. Sec. IV. Although the decay
of total magnetization seems to be slowing down, the data is not conclusive to
determine its long-time behavior. We notice the oscillatory q.s. regime for
$\alpha\lesssim 4$ in Fig. 14 for the single-site observable near the edge,
which was also observed for nonintegrable n.n.n. TFIC Dağ and Sun (2021).
Although the oscillations grow as we decrease $\alpha$, one could still
observe the onset of a q.s. regime in the observable near the edge of the
chain for $\alpha=3$, compared to the observable in the middle of the chain
(Fig. 14b). The difference in the single-site magnetization of spins near the
edge and in the middle of the chain decreases as we keep decreasing $\alpha$.
Fig. 14c, for $\alpha=2.5$, demonstrates for the first time an onset of a q.s.
regime for a spin in the middle of the chain. However, the onset of the q.s.
regime is delayed in the middle of the chain compared to near the edge, and
this observation points to the locality of the underlying Hamiltonian that is
still preserved to an extent. In contrast, there is no boundary effect
observed at $\alpha=2$ (Fig. 14d), where we do not see any difference in the
general trend of the nonequilibrium responses. This is physically intuitive,
because the model possesses algebraic light cones for $\alpha>2$ instead of a
logarithmic cone, pointing to the importance of the locality in the formation
of a q.s. regime. As discussed previously, the origin of the q.s. regime could
be traced back to the asymmetric motion of quasi-particles in opposite
directions, and the reflection from the nearest boundary. When the model is no
longer integrable, even though quasi-particles do not exist one could still
refer to correlation speeds measured through linear or sublinear light cones.
In conclusion, the presence of a q.s. regime does not depend on the
integrability, but rather depends on the locality of the underlying
Hamiltonian. In this sense, this boundary induced temporal regime acts as a
signature of linear, nearly-linear or sublinear light cones.

(a)

(b)

Figure 15: (a) Single-site magnetization at $r=6$ for different
$\alpha=2,4,10$ and different system sizes $N=32,64$ at transverse field
$h=0.5\tilde{h}_{c}$. (b) Total magnetization $\sum_{1}^{N}\sigma_{r}^{z}/N$,
quenched from a polarized state to
$h/\tilde{h}_{c}=1.04,1.05,1.06,1.07,1.08,1.09,1.1$, from top to bottom. The
dashed line-style indicates curves which cross the time axis. The last
magnetization that does not cross the time axis is $h/\tilde{h}_{c}=1.06$,
meaning that the DCP should lie in the interval of $(1.06,1.07)\tilde{h}_{c}$.

To demonstrate that the q.s. regime is not a finite-size effect, we show in
Fig. 15a the single-site magnetizations of a spin close to the boundary at
$r=6$ for different system sizes and different $\alpha$. For a given $\alpha$,
one can determine the time at which finite-size effects kick in by observing
when the data for different system sizes $N=32,64$ no longer overlap. Note
that for $\alpha=10$ and $\alpha=4$, the q.s. regime develops before the
finite-size effects appear.

Finally, we set $\alpha=2.5$, which yields a long-range model Dutta and
Bhattacharjee (2001) with algebraic light cones Foss-Feig  _et al._ (2015),
and study the resulting local order profiles due to the q.s. regime. In all of
our $t$-DMRG calculations on long-range TFIC, we Kac-normalize the interaction
term by dividing it by $\sum_{1}^{L}1/r^{\alpha}$, and hence approximate the
QCP as $\tilde{h}_{c}=J$, which in the thermodynamic limit is a good
approximation when $\alpha>1$ Jaschke  _et al._ (2017). This is performed
because QCP is actually not well-defined for a finite-size system, and finite-
size fluctuations are expected to shift the “critical point” to smaller
values. Therefore, we first determine the DCP based on the total magnetization
Halimeh  _et al._ (2021). Fig. 15b shows the total magnetization for various
transverse field strengths. One sees that the total magnetization crosses the
time axis within the interval $h_{dc}/\tilde{h}_{c}\in(1.06,1.07)$, which pins
down the DCP of the model quenched from a polarized state Halimeh  _et al._
(2021). The observation that $h_{dc}>\tilde{h}_{c}$ points to the fact that
the QCP is initially underestimated.

(a)

(b)

Figure 16: (a) Single-site order profiles calculated with $t$-DMRG for the
long-range TFIC, with $\alpha=2.5$ and system size $N=32$, at different sites
$r=1,3,6$ when a UV temporal cutoff of $t^{*}=2$ and an infrared temporal
cutoff of $t_{l}=6$ are applied. (b) Finite-time analysis curves for $r=3$ are
calculated with $t^{*}=0$ and $t_{l}=2,3,4,5,6$, from top to bottom. Since
there is no crossing point, finite-time analysis fails.

Since the accessible simulation times are significantly limited, we set the
infrared and ultraviolet cutoffs as $t_{l}=6$ and $t^{*}=2$ (Fig. 16a).
Remarkably, we find that the single-site order profiles at different sites dip
at the same $h/\tilde{h}_{c}$ value, suggesting a crossover in the interval of
$h_{cr}/\tilde{h}_{c}\in(1.15,1.2)$. We note that unlike in previous sections
and Ref. Dağ  _et al._ (2021) where we observe
$\overline{|C_{1}(t)|}>\overline{|C_{3}(t)|}>\overline{|C_{6}(t)|}>...$ in the
ordered phase, here for the long-range model we do not observe such an
ordering. This is because of the presence of algebraic, instead of linear,
light cones for $\alpha=2.5$, which wash away the differences in the single-
site magnetization of nearby sites, and hence make the system less locally
connected. Strengthening this argument, we also find that the finite-time
analysis fails for this model. Fig. 16b shows nonequilibrium responses at
$r=3$ between a UV cutoff $t^{*}=0$ and various infrared cutoffs
$t_{l}=2,3,4,5,6$, from top to bottom. The curves do not intersect at any
point, unlike what we observed for n.n.n. TFIC models.

(a)

(b)

(c)

Figure 17: Order profiles calculated with $t$-DMRG for the long-range TFIC at
different system sizes $N=16,24,32$ with $\alpha=2.5$ (a) at site $r=1$ and
$h/\tilde{h}_{c}=0.5$ with respect to different infrared cutoffs parametrized
as $t_{l}(N)/R_{N}$; (b) at $r=1$ with respect to $h/\tilde{h}_{c}$ with an
infrared cutoff of $t_{l}(N,h)$; (c) at $r=3$ with respect to
$h/\tilde{h}_{c}$ with an infrared cutoff of $t_{l}(N,h)$. All subfigures use
a UV cutoff of $t^{*}=0$.

The crossover region might be probing the actual QCP $h_{c}$, as it appears
above the DCP determined with total magnetization at
$h_{dc}/\tilde{h}_{c}\in(1.06,1.07)$, and this is expected in long-range
models for quenches from the ferromagnetic phase Halimeh and Zauner-Stauber
(2017); Homrighausen  _et al._ (2017); Lang  _et al._ (2018). To test this
idea, we apply finite-size analysis with $N=16$ and $N=24$. This analysis
accentuates the shortness of accessible times as the system size increases: A
system size-dependent infrared cutoff $t_{l}(N)=6R_{N}$, where $R_{N}=N/32$,
corresponds to $t_{l}=4.5$ and $t_{l}=3$ for $N=24$ and $N=16$, respectively,
throughout the entire set of transverse-field values. An additional challenge
in the long-range models is the oscillations with large amplitudes Halimeh
_et al._ (2017) as seen in Fig. 14c. This feature particularly becomes
significant as we approach the boundary. Fig. 17a demonstrates how the order
profile at $h/\tilde{h}_{c}=0.5$ would change with different infrared cutoffs
$t_{l}(N)/R_{N}$ and system sizes. The order profiles are also affected by the
oscillations, in particular for early infrared cutoffs, e.g.,
$t_{l}(N)=6R_{N}$. As the averaging interval increases with increasing
$t_{l}$, we notice that the oscillations are damped and the dependence on the
cutoff becomes less pronounced. Again due to oscillations, comparison between
the order profiles of different system sizes becomes less reliable, as can
also be seen in Fig. 17a where the ordering between different system sizes
keeps changing as we change $t_{l}(N)/R_{N}$. While the trend of the data
suggests that the oscillations continue to be damped and the ordering possibly
converges as we increase the infrared cutoff, we cannot confirm this
conclusion due to limited simulation times. Therefore, both to alleviate the
effects of oscillations and to utilize the accessible times as much as
possible, we adopt an infrared cutoff that is a function of both system size
and the transverse field $t_{l}(N,h)=R_{N}t_{max}(N=32,h)$ where
$t_{max}(N=32,h)$ is the longest time accessible in our $t$-DMRG calculations
for the largest considered system size $N=32$ at each $h$. While this adaptive
approach helps to decrease the effect of oscillations on the order profiles,
in particular deep in the ferromagnetic phase, it does not completely
eliminate the effect, which is possible only by increasing the simulation
times. Figs. 17b and 17c show the order profiles of different sizes at $r=1$
and $r=3$, respectively. We observe that the order profile of $N=32$ tends to
decrease at much smaller $h$ than where the crossover resides. This is exactly
because of the limited simulation times and hence the application of very
early infrared cutoffs as we increase the transverse field. It is also worth
noting that finite-size scaling for models with power-law interaction profiles
faces the additional challenge that the profile tails will look different for
different system sizes, and this becomes particularly nontrivial near
criticality. As a result of the above, whether the observed crossover profile
is actually a transition cannot be answered conclusively.

##  VI Conclusions and Outlook

We have studied the q.s. regime at sites close to the edge in both integrable
and nonintegrable TFICs, and shown the necessity of the zero modes to
stabilize this temporal regime. Based on analytical arguments and numerical
evidence, we demonstrated the presence of non-analyticity in the dynamical
order profiles of integrable TFIC and near-integrable TFIC treated with qMFT.
This suggested a dynamical critical point that coincides with the QCP in the
integrable TFIC and slightly shifts from the corresponding QCP in the near-
integrable model. Therefore, we conclude that the single-site observables are
able to extract QCP in short-range interacting models independently of
integrability, measurement location and initial state so long as $r\ll N/2$
and $h_{i}<h_{c}$.

We performed $t$-DMRG calculations to obtain the order profiles of single-site
observables at various sites in open-boundary nonintegrable TFICs. The models
that we studied range from locally connected with next nearest neighbor
interactions to power-law decaying interactions. Our qMFT and $t-$DMRG results
for the location of the DCP in near-integrable model agree, as we found
$h_{dc}^{\Delta=0.1}\in(1.1,1.15)$ with $t-$DMRG for system sizes $N\leq 120$.
Finite-size analysis for the strongly nonintegrable model determined the DCP
to be in the interval of $h_{dc}^{\Delta=1}\in(2.3,2.35)$, independent of the
measurement sites, and to be shifted from the QCP at $h_{c}\approx 2.46$.
Whether this shift is related to finite-size effects, or due to interactions
as argued in Titum  _et al._ (2019), is a question that needs further
elaboration in the future. We find that the finite-time analysis with single-
site observables matches well with the finite-size analysis in the near-
integrable model, whereas there is a mismatch between the two analyses in the
strongly nonintegrable model. We argue that this is likely because of a change
in the light cone structure as the transverse field increases, which is yet to
be explored. We demonstrated the presence of a q.s. regime even in long-range
power-law interacting TFIC for $\alpha=3$ and $\alpha=2.5$, where in the
former one could still observe a significant difference in the single-site
magnetization of the spin in the middle and close to an edge, pointing to
nearly-linear light cones hosted in the system. Increasing the range of
interactions to $\alpha=2.5$ decreases the differences between sites,
emphasizing the long-range nature of the model and expectantly the finite-time
analysis fails. The order profiles at $\alpha=2.5$ for different sites dip at
the same transverse field $h_{cr}/\tilde{h}_{c}<1.2$ suggesting a crossover,
which is found to be larger than the calculated DCP based on total
magnetization at $h_{dc}/\tilde{h}_{c}=(1.06,1.07)$. Given that the crossover
dip appears after $h_{dc}/\tilde{h}_{c}$, and since the calculations were
based on an estimate of QCP $\tilde{h}_{c}$, it is an interesting direction to
check whether the crossover point $h_{cr}$ probes the actual QPT. Our data is
significantly limited to early times, especially as we approach the transition
boundary, and hence is currently inconclusive to answer this question.

Finally, our setup is experimentally convenient, because (i) open-boundary
chains are a more natural setup than their periodic counterparts in most
quantum simulators, and (ii) spatially minimal probes are readily accessible
in modern quantum simulators Bakr  _et al._ (2009). Most theoretical works
have naturally focused on infinite or periodic chains to utilize the
translational symmetry Halimeh  _et al._ (2017); Titum  _et al._ (2019);
Halimeh  _et al._ (2020), which removes site-dependency of the local order
profiles within the ordered phase. In this sense, our work complements the
literature via explicitly demonstrating the potential of spatially minimal
measurements in open-boundary chains and hence exploiting the boundary effects
in probing criticality.

###### Acknowledgements.

The authors thank Philipp Hauke for suggestions on the manuscript. C.B.D.
thanks Kai Sun for extensive and fruitful discussions; Norman Yao and L.-M.
Duan for helpful suggestions. C.B.D. was supported by the National Science
Foundation under Grant EFRI-1741618 and ITAMP grant at Harvard University.
This work is part of and supported by Provincia Autonoma di Trento, the ERC
Starting Grant StrEnQTh (project ID 804305), the Google Research Scholar Award
ProGauge, and Q@TN — Quantum Science and Technology in Trento. I.P.M.
acknowledges support from the Australian Research Council (ARC) Discovery
Project Grants No. DP190101515 and DP200103760.

## References

  * Landau (1937) L. D. Landau, Zh. Eksp. Teor. Fiz. 11, 19 (1937).
  * Cardy  _et al._ (1996) J. Cardy, P. Goddard, and J. Yeomans,  _Scaling and Renormalization in Statistical Physics_, Cambridge Lecture Notes in Physics (Cambridge University Press, 1996).
  * Kadanoff (1993) L. P. Kadanoff, “Scaling and universality in statistical physics,” in  _From Order to Chaos_ (World Scientific, 1993) pp. 280–293.
  * Zinn-Justin (2002) J. Zinn-Justin,  _Quantum Field Theory and Critical Phenomena_, Internat. Ser. Mono. Phys. (Clarendon Press, Oxford, 2002).
  * Sachdev (2001) S. Sachdev,  _Quantum Phase Transitions_ (Cambridge University Press, 2001).
  * Hohenberg and Halperin (1977) P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Mori  _et al._ (2018) T. Mori, T. N. Ikeda, E. Kaminishi, and M. Ueda, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 112001 (2018).
  * Heyl (2018) M. Heyl, Reports on Progress in Physics 81, 054001 (2018).
  * Calabrese and Cardy (2006) P. Calabrese and J. Cardy, Phys. Rev. Lett. 96, 136801 (2006).
  * Calabrese  _et al._ (2011) P. Calabrese, F. H. L. Essler, and M. Fagotti, Phys. Rev. Lett. 106, 227203 (2011).
  * Kollath  _et al._ (2007) C. Kollath, A. M. Läuchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007).
  * Calabrese  _et al._ (2012) P. Calabrese, F. H. L. Essler, and M. Fagotti, Journal of Statistical Mechanics: Theory and Experiment 2012, P07016 (2012).
  * Tsuji  _et al._ (2013) N. Tsuji, M. Eckstein, and P. Werner, Phys. Rev. Lett. 110, 136404 (2013).
  * Heyl  _et al._ (2013) M. Heyl, A. Polkovnikov, and S. Kehrein, Phys. Rev. Lett. 110, 135704 (2013).
  * Sciolla and Biroli (2013) B. Sciolla and G. Biroli, Phys. Rev. B 88, 201110 (2013).
  * Chiocchetta  _et al._ (2015) A. Chiocchetta, M. Tavora, A. Gambassi, and A. Mitra, Phys. Rev. B 91, 220302 (2015).
  * Halimeh and Zauner-Stauber (2017) J. C. Halimeh and V. Zauner-Stauber, Phys. Rev. B 96, 134427 (2017).
  * Karl  _et al._ (2017) M. Karl, H. Cakir, J. C. Halimeh, M. K. Oberthaler, M. Kastner, and T. Gasenzer, Phys. Rev. E 96, 022110 (2017).
  * Titum  _et al._ (2019) P. Titum, J. T. Iosue, J. R. Garrison, A. V. Gorshkov, and Z.-X. Gong, Phys. Rev. Lett. 123, 115701 (2019).
  * Gurarie (2019) V. Gurarie, Phys. Rev. A 100, 031601 (2019).
  * Uhrich  _et al._ (2020) P. Uhrich, N. Defenu, R. Jafari, and J. C. Halimeh, Phys. Rev. B 101, 245148 (2020).
  * Wu (2020a) Y. Wu, Phys. Rev. B 101, 014305 (2020a).
  * Wu (2020b) Y. Wu, Phys. Rev. B 101, 064427 (2020b).
  * Dağ and Sun (2021) C. B. Dağ and K. Sun, Phys. Rev. B 103, 214402 (2021).
  * Halimeh  _et al._ (2021) J. C. Halimeh, D. Trapin, M. Van Damme, and M. Heyl, Phys. Rev. B 104, 075130 (2021).
  * Haldar  _et al._ (2021) A. Haldar, K. Mallayya, M. Heyl, F. Pollmann, M. Rigol, and A. Das, Phys. Rev. X 11, 031062 (2021).
  * Halimeh  _et al._ (2019) J. C. Halimeh, N. Yegovtsev, and V. Gurarie, (2019), arXiv:1903.03109 [cond-mat.stat-mech] .
  * Trapin  _et al._ (2021) D. Trapin, J. C. Halimeh, and M. Heyl, Phys. Rev. B 104, 115159 (2021).
  * Heyl  _et al._ (2018) M. Heyl, F. Pollmann, and B. Dóra, Phys. Rev. Lett. 121, 016801 (2018).
  * Dağ  _et al._ (2019) C. B. Dağ, K. Sun, and L.-M. Duan, Phys. Rev. Lett. 123, 140602 (2019).
  * Gómez-Ruiz  _et al._ (2018) F. J. Gómez-Ruiz, J. J. Mendoza-Arenas, F. J. Rodríguez, C. Tejedor, and L. Quiroga, Phys. Rev. B 97, 235134 (2018).
  * Dağ  _et al._ (2020) C. B. Dağ, L.-M. Duan, and K. Sun, Phys. Rev. B 101, 104415 (2020).
  * Sun  _et al._ (2020) Z.-H. Sun, J.-Q. Cai, Q.-C. Tang, Y. Hu, and H. Fan, Annalen der Physik 532, 1900270 (2020).
  * Bandyopadhyay  _et al._ (2021a) S. Bandyopadhyay, A. Polkovnikov, and A. Dutta, Phys. Rev. Lett. 126, 200602 (2021a).
  * Bandyopadhyay  _et al._ (2021b) S. Bandyopadhyay, P. Uhrich, A. Paviglianiti, and P. Hauke, (2021b), arXiv:2108.01718 [cond-mat.str-el] .
  * Nicklas  _et al._ (2015) E. Nicklas, M. Karl, M. Höfer, A. Johnson, W. Muessel, H. Strobel, J. Tomkovič, T. Gasenzer, and M. K. Oberthaler, Phys. Rev. Lett. 115, 245301 (2015).
  * 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, 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, Nature 551, 601–604 (2017).
  * Hauke  _et al._ (2014) P. Hauke, M. Lewenstein, and A. Eckardt, Phys. Rev. Lett. 113, 045303 (2014).
  * 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, Nature Physics 14, 265 (2018).
  * Yang  _et al._ (2019) H.-X. Yang, T. Tian, Y.-B. Yang, L.-Y. Qiu, H.-Y. Liang, A.-J. Chu, C. B. Dağ, Y. Xu, Y. Liu, and L.-M. Duan, Phys. Rev. A 100, 013622 (2019).
  * Muniz  _et al._ (2020) J. A. Muniz, D. Barberena, R. J. Lewis-Swan, D. J. Young, J. R. K. Cline, A. M. Rey, and J. K. Thompson, Nature 580, 602 (2020).
  * Nie  _et al._ (2020) X. Nie, B.-B. Wei, X. Chen, Z. Zhang, X. Zhao, C. Qiu, Y. Tian, Y. Ji, T. Xin, D. Lu, and J. Li, Phys. Rev. Lett. 124, 250601 (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, Science Advances 6 (2020), 10.1126/sciadv.aba4935.
  * Chen  _et al._ (2020) B. Chen, X. Hou, F. Zhou, P. Qian, H. Shen, and N. Xu, Applied Physics Letters 116, 194002 (2020), https://doi.org/10.1063/5.0004152 .
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Bakr  _et al._ (2009) W. S. Bakr, J. I. Gillen, A. Peng, S. Fölling, and M. Greiner, Nature (London) 462, 74 (2009), arXiv:0908.0174 [cond-mat.quant-gas] .
  * Dağ  _et al._ (2018) C. B. Dağ, S.-T. Wang, and L.-M. Duan, Phys. Rev. A 97, 023603 (2018).
  * Dağ  _et al._ (2021) C. B. Dağ, Y. Wang, P. Uhrich, X. Na, and J. C. Halimeh, (2021), arXiv:2105.05986 [cond-mat.quant-gas] .
  * Iglói and Rieger (2011) F. Iglói and H. Rieger, Phys. Rev. Lett. 106, 035701 (2011).
  * Foss-Feig  _et al._ (2015) M. Foss-Feig, Z.-X. Gong, C. W. Clark, and A. V. Gorshkov, Phys. Rev. Lett. 114, 157201 (2015).
  * Hauke and Tagliacozzo (2013) P. Hauke and L. Tagliacozzo, Phys. Rev. Lett. 111, 207202 (2013).
  * Budich and Heyl (2016) J. C. Budich and M. Heyl, Phys. Rev. B 93, 085416 (2016).
  * Lieb  _et al._ (1961) E. Lieb, T. Schultz, and D. Mattis, Annals of Physics 16, 407 (1961).
  * Wang  _et al._ (2021) Z. Wang, M. Foss-Feig, and K. R. A. Hazzard, Physical Review Research 3 (2021), 10.1103/physrevresearch.3.l032047.
  * Yang (1952) C. N. Yang, Phys. Rev. 85, 808 (1952).
  * Iglói  _et al._ (2013) F. Iglói, G. Roósz, and Y.-C. Lin, New Journal of Physics 15, 023036 (2013).
  * Gogolin and Eisert (2016) C. Gogolin and J. Eisert, Reports on Progress in Physics 79, 056001 (2016).
  * Kemp  _et al._ (2017) J. Kemp, N. Y. Yao, C. R. Laumann, and P. Fendley, Journal of Statistical Mechanics: Theory and Experiment 2017, 063105 (2017).
  * (61) ITensor Library (version 2.0.11) http://itensor.org .
  * Neyenhuis  _et al._ (2017) B. Neyenhuis, J. Zhang, P. W. Hess, J. Smith, A. C. Lee, P. Richerme, Z.-X. Gong, A. V. Gorshkov, and C. Monroe, Science Advances 3 (2017), 10.1126/sciadv.1700672.
  * Dutta and Bhattacharjee (2001) A. Dutta and J. K. Bhattacharjee, Phys. Rev. B 64, 184106 (2001).
  * Knap  _et al._ (2013) M. Knap, A. Kantian, T. Giamarchi, I. Bloch, M. D. Lukin, and E. Demler, Phys. Rev. Lett. 111, 147205 (2013).
  * Schollwöck (2005) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
  * Schollwöck (2011) U. Schollwöck, Annals of Physics 326, 96 (2011), january 2011 Special Issue.
  * Daley  _et al._ (2004) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Journal of Statistical Mechanics: Theory and Experiment 2004, P04005 (2004).
  * Gobert  _et al._ (2005) D. Gobert, C. Kollath, U. Schollwöck, and G. Schütz, Phys. Rev. E 71, 036102 (2005).
  * (69) Ian P. McCulloch, Matrix Product Toolkit. URL: http://physics.uq.edu.au/people/ianmcc/mptoolkit/.
  * (70) A. N. Krylov, “On the numerical solution of the equation by which in technical questions frequencies of small oscillations of material systems are determined,” Izvestija AN SSSR (News of Academy of Sciences of the USSR), Otdel. mat. i estest. nauk VII, 491–539 (1931), in Russian.
  * Crosswhite  _et al._ (2008) G. M. Crosswhite, A. C. Doherty, and G. Vidal, Phys. Rev. B 78, 035116 (2008).
  * Jaschke  _et al._ (2017) D. Jaschke, K. Maeda, J. D. Whalen, M. L. Wall, and L. D. Carr, New Journal of Physics 19, 033032 (2017).
  * Homrighausen  _et al._ (2017) I. Homrighausen, N. O. Abeling, V. Zauner-Stauber, and J. C. Halimeh, Phys. Rev. B 96, 104436 (2017).
  * Lang  _et al._ (2018) J. Lang, B. Frank, and J. C. Halimeh, Phys. Rev. B 97, 174401 (2018).
  * Halimeh  _et al._ (2017) J. C. Halimeh, V. Zauner-Stauber, I. P. McCulloch, I. de Vega, U. Schollwöck, and M. Kastner, Phys. Rev. B 95, 024302 (2017).
  * Halimeh  _et al._ (2020) J. C. Halimeh, M. Van Damme, V. Zauner-Stauber, and L. Vanderstraeten, Phys. Rev. Research 2, 033111 (2020).
  * Vekić and White (1993) M. Vekić and S. R. White, Phys. Rev. Lett. 71, 4283 (1993).

##  Appendix A Independency of the results from temporal cutoffs

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

Figure 18: The single-site nonequilibrium phase diagrams of the integrable
TFIC with a UV temporal cutoff of $t^{*}=10$ for (a) $C_{3}(t)$, (b)
$C_{9}(t)$ and (c) $C_{12}(t)$; with a UV temporal cutoff of $t^{*}=20$ for
(d) $C_{3}(t)$, (e) $C_{6}(t)$, (f) $C_{9}(t)$ and (g) $C_{12}(t)$; with a UV
temporal cutoff of $t^{*}=2\alpha\Delta x/v_{q}$ where $\Delta x=r-1$ with $r$
being the single-site observable location and $\alpha=2$ is a tuning parameter
for (h) $C_{3}(t)$, (i) $C_{6}(t)$, (j) $C_{9}(t)$ and (k) $C_{12}(t)$. (l-o)
The single-site nonequilibrium phase diagrams with an initial state as the
ground state of an initial Hamiltonian with $h_{i}=0.1$ and a UV temporal
cutoff of $t^{*}=10$ for (l) $C_{3}(t)$, (m) $C_{6}(t)$, (n) $C_{9}(t)$ and
(o) $C_{12}(t)$. The behavior is qualitatively the same as the results of
$h_{i}=0$.

We first test whether our results depend on the choice of ultraviolet (UV)
cutoff $t^{*}$. The results in the main text are produced with a fixed UV
cutoff of $t^{*}=10$, for all $h$. However, none of our results depend on the
choice of UV cutoff: Subfigures (a-k) in Fig. 18 all show the same qualitative
behavior for single-site nonequilibrium phase diagrams for various choices of
UV cutoff. Figs. 18a – 18c complement the $t^{*}=10$ data in the main text by
showing the nonequilibrium phase diagrams of the single-site magnetization at
$r=3,9,12$. Figs. 18d – 18g exhibit another fixed UV cutoff of $t^{*}=20$,
whereas Figs. 18h – 18k demonstrate the results of a parametric UV cutoff for
all studied sites. The parametric UV cutoff is determined as follows: We
roughly estimate the onset of the q.s. regime as the time required for the
quasiparticles to reflect back from the edge closest to the observation site.
Therefore, the estimate can be mathematically stated as $t^{*}=2\alpha\Delta
x/v_{q}$, where the distance $\Delta x=r-1$ is the distance between the
observation site, $r=3,6,9,12$ and the closest edge site, $r^{\prime}=1$, in
our case. The parameter $\alpha$ is a tuning parameter, as our analytical
formula is only an estimate. In fact we find that $\alpha=2$ presents phase
diagrams qualitatively the same as those of other $t^{*}$, for all $r$ that we
studied.

The long-time asymptotic value of $J_{1}(4t)/2t$ is $t^{-3/2}$, which is what
we observe in Fig. 5a not only for the edge but for all probe sites $r$. Let
us note that the same holds for the near-integrable model treated with qMFT
(Fig. 8c) as well as the disordered phase in both models (Figs. 5b and 8d).
Therefore, it is straightforward to find how the time-averaged magnetization,
e.g., single-site dynamical order parameters, scale with the system size at
the QCP and in the disordered phase.

$\displaystyle\overline{C(t,h\geq h_{c})}$ $\displaystyle=$
$\displaystyle\frac{1}{(N/v_{c}-t^{*})}\int_{t^{*}}^{N/v_{c}}t^{-3/2}dt=-2\frac{1}{(N/v_{c}-t^{*})}t^{-1/2}\bigg{|}_{t^{*}}^{N/v_{c}}=-2\frac{1}{(N/v_{c}-t^{*})}\left(\frac{1}{\sqrt{N/v_{c}}}-\frac{1}{\sqrt{t^{*}}}\right)$
(62) $\displaystyle=$
$\displaystyle-2\frac{(4h_{c})^{3/2}}{N^{3/2}-4t^{*}h_{c}\sqrt{N}}+2\frac{4h_{c}}{(N-4t^{*}h_{c})\sqrt{t^{*}}}.$

Note that $t^{*}$, the ultraviolet cutoff, has to be introduced in the first
line to prevent divergence which is nonphysical. Series expansion for Eq. (62)
when $N\rightarrow\infty$ gives,

$\displaystyle\overline{C(t,h\geq
h_{c})}=\frac{8}{\sqrt{t^{*}}}\left(\frac{N}{h_{c}}\right)^{-1}-16\left(\frac{N}{h_{c}}\right)^{-3/2}+32\sqrt{t^{*}}\left(\frac{N}{h_{c}}\right)^{-2}+\mathcal{O}(N^{-5/2}).$
(63)

According to Eq. (63), the dynamical order at the QCP decreases as $\propto
1/N$ for all probe sites $r$ and becomes $0$ in the thermodynamic limit. We
numerically confirm this prediction in Figs. 19a, 19b and 19c with fits of the
form $N^{-b}$ where $b\approx 1$. Let us note that the analytical expression
is only exact asymptotically in time and the numerics take into account the
early-time behavior of $|C_{r}(t)|$ too. Hence, it is not surprising to see
$b\approx 1$ instead of $b=1$ in the numerics.

(a)

(b)

(c)

(d)

Figure 19: The system size scaling of different single-site magnetization at
the critical point $h_{c}$ with an initial state at $h_{i}=0$ and cutoff (a)
$t^{*}=10$, (b) $t^{*}=20$, (c) $t^{*}=2\alpha\Delta x/v_{q}$ and $\alpha=2$,
and (d) for an initial state at $h_{i}=0.1$ and a UV temporal cutoff of
$t^{*}=10$. A downward trend can be seen with $aN^{-b}$ where $b\approx 1$ in
all curves. In each subfigure, the curves from up to bottom represent the
single sites $r=3,6,9,12$, respectively, as seen in the legend.

(a)

(b)

(c)

(d)

Figure 20: The energy gap analysis (a-b) for the near integrable model with
DMRG, (c-d) for the effective noninteracting Hamiltonian derived through qMFT
method. (a,c) The critical points are marked $h_{c}=1.1631\pm 0.0037$ and
$h_{c}^{\text{qmft}}=1.1438\pm 0.0011$, respectively. (b,d) The energy gap
scales with system size as $1/N$ in both.

##  Appendix B Initial state independence of the results

In this subsection, we change the initial state to the ground state of an
initial Hamiltonian with $h_{i}=0.1$, and test whether any of our results
depend on the initial state. Figs. 18l – 18o show the single-site
nonequilibrium phase diagrams computed with this initial state. We do not
observe a change in the qualitative behavior. The single-site magnetization at
the critical point still decreases with increasing system size in Fig. 19d by
exhibiting $b\approx 1$.

The analytical expression for the q.s. value of the edge magnetization in the
ordered phase is,

$\displaystyle
C^{qs}_{1}(h,h_{i})=\frac{(1-h^{2})(1-2h_{i})^{1/2}}{1-2hh_{i}},$ (64)

for $h,h_{i}<1$. Let us rewrite it in terms of the reduced control parameter
$h_{n}=(h_{c}-h)/h_{c}$ as

$\displaystyle
C^{qs}_{1}(h_{n},h_{i})=\frac{(2-h_{n})h_{n}(1-2h_{i})^{1/2}}{1+2(h_{n}-1)h_{i}}.$
(65)

In the vicinity of the transition, $h_{n}\rightarrow 0$, we can expand this
expression and find

$\displaystyle C^{qs}_{1}(h_{n}\rightarrow 0,h_{i})$ $\displaystyle=$
$\displaystyle a_{1}(h_{i})h_{n}+\mathcal{O}(h_{n}^{2})$ (66) $\displaystyle
a_{1}(h_{i})$ $\displaystyle=$ $\displaystyle
2+h_{i}+\frac{3}{4}h_{i}^{2}+\frac{5}{8}h_{i}^{3}+\cdots,$

Therefore, one can see that $C^{qs}_{1}(h_{n})\propto h_{n}$ in the vicinity
of the transition, $h_{n}\rightarrow 0$, regardless of the choice of initial
state. The initial state only changes the coefficient in front of $h_{n}$,
which is nonuniversal.

##  Appendix C Energy gap analysis

Here we present the ground state energy gap analysis of the near integrable
model $\Delta=0.1$ calculated with DMRG. Then we discuss the single-particle
energy gap analysis for the effective Hamiltonian produced by the qMFT method.

Figs. 20a shows how the location of the ground state energy gap minimum scales
with system size. This finite-size scaling analysis computed with DMRG gives
$h_{c}\approx 1.1631$. Meanwhile the energy gap at the QCP scales as $N^{-1}$
which is shown in Fig. 20b.

qMFT method provides us an effective Hamiltonian for the near-integrable
model. In Fig. 20c we focus on the minimum of the energy gap of this
Hamiltonian. We observe that location of the minimum scales with the system
size giving $h_{c}^{\text{qmft}}\approx 1.14378$. The energy gap at
$h_{c}^{\text{qmft}}$ scales with the system size as $N^{-1}$. This result is
consistent with the observation that a critically prethermal regime appears as
we suddenly quench to the vicinity of the DCP of the near-integrable model in
Ref. Dağ  _et al._ (2021). Because the DCP of the near-integrable model,
$h_{dc}\approx 1.1437$ is almost equal to the QCP of the effective qMFT
Hamiltonian. Then the duration of the prethermal regime is actually governed
by the QCP dynamics of the effective qMFT Hamiltonian in Dağ  _et al._ (2021).

##  Appendix D Changing the boundary conditions

Figure 21: Different boundary conditions are compared for the integrable TFIC
at $h=0.5$ and system size $N=24$. One can obtain a q.s. regime with smooth
boundary conditions too. OBC, PBC and SBC stand for open, periodic and smooth
boundary conditions, respectively.

We demonstrate that the q.s. temporal regime emerges not only when we
introduce hard boundaries Iglói and Rieger (2011), but also for smooth
boundaries. A smooth boundary condition can be applied by smoothly turning off
the Hamiltonian parameters towards the edges of the chain Vekić and White
(1993). Fig. 21 shows the single-site nonequilibrium responses of the
integrable TFIC with hard boundarieas (red-diamonds), smooth boundaries
(green-solid) and periodic boundary condition (blue-dotted). As shown before
Dağ and Sun (2021), the middle of a hard-boundary chain (yellow triangles)
acts like an arbitrary site in a periodic chain. We thus conclude that the
q.s. regime is robust against altering the boundary conditions, so long as
they remain open.

##  Appendix E Benchmarking the qMFT method

(a)

(b)

(c)

Figure 22: Benchmarking qMFT analysis. All subfigures compare the results of
qMFT, $t$-DMRG and exact diagonalization (ED) results (see individual legends
for system sizes) for $C_{3}(t)$. The external fields are (a) $h=0.5$, (b)
$h=1.1$ and (c) $h=1.2$. In all subfigures, the qMFT nonequilibrium response
matches sufficiently well with the nonequilibrium responses of the exact
methods, in particular in early times close to QPT and in the disordered phase
where exact methods are constrained either by finite-size (ED algorithm) or
finite bond dimension and short simulation times (t-DMRG algorithm).

(a)

(b)

Figure 23: Non-equilibrium responses calculated with t-DMRG for weakly
interacting TFIM with $\Delta/J=0.1$ at $N=96$ are compared with respect to
different maximum bond dimensions, $\chi_{m}=100$ and $\chi_{m}=200$ which is
introduced as an approximation in the t-DRMG algorithm to simulate for longer
times, see Sec. IV.3. The comparison is shown in subfigures (a) at $h/J=1.3$
and (b) at $h/J=1$. (a) With increasing bond dimension, a downward moving
envelope appears in the t-DMRG result, emphasized by the dotted-gray arrow,
agreeing with the qMFT results in the disordered phase. The solid-black arrow
points to the effects of finite bond dimension. (b) The qMFT result at
$N=1440$ is compared with the t-DMRG results. We notice the convergence of the
t-DMRG result with greater $\chi_{m}$ (red-solid) to the qMFT result (green-
dotted) for a small time interval that was accessible, before showing the
effects of finite bond dimension.

Fig. 22 compares the results of qMFT, $t$-DMRG and exact diagonalization (ED)
at $h=0.5$ and $h=1.1$ in the ordered phase and $h=1.2$ in the disordered
phase. We observe that the qMFT analysis can even capture the correct
frequency of the oscillations in early times and the general trend of the
nonequilibrium response successfully, although it does not totally match with
the exact methods, which is expected due to the fact that it is an approximate
method that averages out the interactions. In Fig. 22a, both the ED and t-DMRG
responses show finite-size effects in the form of quantum revivals (solid-blue
and red-dashed-dotted). The finite-size effects of these methods can also be
clearly seen in Figs. 22b and 22c. As we approach the QCP and in the
disordered phase, the match between the exact methods and the qMFT does not
survive past $t=10$. This is either because of the finite-size effects, as
apparent in ED response, or the finite maximum bond dimensions set in the
t-DMRG algorithm, see Sec. IV.3 for more details. In order to demonstrate the
effects of finite bond dimension, we compare the nonequilibrium responses
calculated with t-DMRG for different maximum bond dimensions. Fig. 23 shows
these responses at $h/J=1.3$ in subfigure (a) and at $h=1$ in subfigure (b).
We observe that the response with bigger maximum bond dimension,
$\chi_{m}=200$ departs from the response with $\chi_{m}=100$ at around
$t\approx 10$, exhibiting a downward envelope and hence agreeing with the qMFT
nonequilibrium responses in the disordered phase, see Fig. 23a. This feature
is shown with a grey-dotted arrow. As time increases, the nonequilibrium
response with $\chi_{m}=200$ starts to exhibit features similar to those of
the response with $\chi_{m}=100$, e.g., recurrences, as pointed out with a
black-solid arrow. We note that these are effects of finite maximum bond
dimension set in our t-DMRG algorithm. One could predict that as the maximum
bond dimension $\chi_{m}$ increases, such effects will occur later in time,
and instead the response will follow closely to that of qMFT response with an
oscillatory downward trend. However, there is a trade-off between the maximum
bond dimension and the maximum accessible simulation time, and simulating
longer times while achieving a satisfactory precision for the response is
simply out of our computational reach. Fig. 23b shows the same comparison
between different $\chi_{m}$ alongside with the qMFT result at $h=1$, which is
a point in the ordered phase. The $t$-DMRG result with $\chi_{m}=100$ shows an
upward trend as time increases, in contrast to the qMFT result which reaches a
plateau, as expected from our theory based on noninteracting fermions. We also
plot the response with $\chi_{m}=200$ and observe that the nonequilibrium
response of greater maximum bond dimension approaches that of the qMFT
results, departing from the t-DMRG result with $\chi_{m}=100$ for a small
interval of time. Additionally we notice that the response with $\chi_{m}=200$
(red-solid) also demonstrates effects of finite maximum bond dimension in
later times with a recurrence similar to what is observed in the response with
smaller bond dimension $\chi_{m}=100$ (blue-star). Therefore, we argue that as
$\chi_{m}$ increases, the t-DMRG result should converge to the qMFT result.
This observation reveals an interesting feature where the approximate qMFT
result can be obtained in full quantum dynamics only with sufficiently large
maximum bond dimensions.