# Spatiotemporal Quenches in Long-Range Hamiltonians

Simon Bernier<EMAIL_ADDRESS>Department of Physics, McGill
University, Montréal, Québec H3A 2T8, Canada  Kartiek Agarwal
<EMAIL_ADDRESS>Department of Physics, McGill University, Montréal,
Québec H3A 2T8, Canada

###### Abstract

Spatiotemporal quenches are efficient at preparing ground states of critical
Hamiltonians that have emergent low-energy descriptions with Lorentz
invariance [1, 2]. The critical transverse field Ising model with nearest
neighbor interactions, for instance, maps to free fermions with a relativistic
low energy dispersion. However, spin models realized in artificial quantum
simulators based on neutral Rydberg atoms, or trapped ions, generically
exhibit long range power-law decay of interactions with $J(r)\sim
1/r^{\alpha}$ for a wide range of $\alpha$. In this work, we study the fate of
spatiotemporal quenches in these models with a fixed velocity $v$ for the
propagation of the quench front, using the numerical time-dependent
variational principle. For $\alpha\gtrsim 3$, where the critical theory is
suggested to have a dynamical critical exponent $z=1$, our simulations show
that optimal cooling is achieved when the front velocity $v$ approaches $c$,
the effective speed of excitations in the critical model. The energy density
is inhomogeneously distributed in space, with prominent hot regions populated
by excitations co-propagating with the quench front, and cold regions
populated by counter-propagating excitations. Lowering $\alpha$ largely blurs
the boundaries between these regions. For $\alpha<3$, we find that the Doppler
cooling effect disappears, as expected from renormalization group results for
the critical model which suggest a dispersion $\omega\sim q^{z}$ with $z<1$.
Instead, we show that excitations are controlled by two relevant length scales
whose ratio is related to that of the front velocity to a threshold velocity
that ultimately determines the adiabaticity of the quench.

††preprint: APS/123-QED

##  I Introduction

Modern quantum simulators hold immense potential for studying fundamental
aspects of quantum many-body systems and materials. Recent experiments in
ultracold atoms [3, 4, 5, 6] and trapped ions [7, 8, 9] have successfully
demonstrated many novel quantum phenomena—a variety of spin models [10, 11,
12, 13, 14, 15], topological quantum numbers [16, 17, 18, 19, 20], many-body
localization [21, 22, 23, 24, 25, 26], lattice gauge theories [27, 28, 29,
30], among others [31, 32, 33]—and have emerged as candidates for programmable
quantum computing [34, 35]. A key application of such artificial quantum
matter is to simulate strongly correlated phases of electrons in conventional
materials [36, 37, 38, 39]. Although Hamiltonians of many such systems can be
approximately realized using a combination of fixed potentials and driving, it
remains a challenge to prepare the system in a state corresponding to a low
enough effective temperature at which the ground state properties can be
reliably explored [36].

Conventionally, state preparation proceeds via adiabatic evolution [40]. The
system is initialized in (or close to) the ground state of a Hamiltonian that
is easy to prepare—the Hamiltonian is usually gapped and the ground state has
low entanglement. The parameters of the Hamiltonian are then tuned such that
the state evolves into the target state, which is often the ground state of a
target Hamiltonian. If the parameters are tuned slowly enough, the quantum
state stays in the ground state of the instantaneous Hamiltonian. However, the
time required to adiabatically prepare a state scales as the square of the
inverse of the smallest energy gap encountered when tuning to the target
Hamiltonian [40]. If the gap closes during evolution, excitations are
inevitably produced and adiabatic techniques fail to produce the target state
with high probability.

In cases where adiabatic evolution fails or takes longer than the coherence
time of the quantum simulator, shortcuts to adiabaticity are required. For
this purpose, counter-diabatic driving was introduced to counter the
production of excitations using auxiliary time-dependent Hamiltonians [41, 42,
43]. Optimal control protocols such as bang-bang protocols have been developed
[44, 45, 46, 35] and rely on classical optimization of the protocol. Spatially
inhomogeneous quenches have also been developed, where portions of a system
act as a sink for excitations [47, 48].

For systems that exhibit emergent Lorentz symmetry, an efficient route to
preparing the ground state of Hamiltonians is via spatiotemporal quenches [49,
50, 1, 2]. This class of protocols can be used to rapidly produce the ground
state of such Hamiltonians even in the critical case, characterized by a
linearly dispersing mode with a minimum energy gap that vanishes as $\sim
1/L$, where $L$ is the linear dimension of the system. In particular, the
system is initialized in a low-entanglement state corresponding to the ground
state of a Hamiltonian that has a gapping perturbation on top of the critical
Hamiltonian. The gapping perturbation is then turned off along a quench front
moving at a time-dependent velocity $v(t)$ greater than the speed of “light”
$c$ of the critical theory. In the simplest version of the protocol, $v(t)$ is
constant in time and optimal cooling is obtained in the limit $v\rightarrow
c^{+}$; see Fig. 1. Such methods should be applicable to quantum simulators
trying to obtain low-energy states of the Hubbard model in two-dimensions [36,
51] (in particular at half filling and large $U$, where a linear spin wave
dispersion emerges) and in one dimensional quantum gases where a low-energy
Luttinger liquid description often applies [52].

Intuitively, the protocol uses Doppler-shifts to result in cooling. In
particular, the quench front excites modes in a chiral way. Modes co-
propagating with the front are blue shifted while counter-propagating modes
are red-shifted. As the velocity of the front approaches the speed of light,
counter-propagating excitations are completely suppressed, and all energy is
carried by excitations propagating along with the quench front, leaving behind
a system with critical ground state correlations. This method thus prepares
the ground state of critical models in a time that scales linearly with system
size, providing a parametric advantage over adiabatic evolution, which
requires a time increasing quadratically with system size [50].

In this work, we study these spatiotemporal quenches with particular emphasis
on an aspect inherent to many artificial simulators based on trapped ion
setups or neutral Rydberg atoms. These systems generically realize effective
spin models with long-range hopping and interactions, with terms decaying as
$1/r^{\alpha}$ in distance $r$, for $0<\alpha\leq 6$. These long-range
interactions (LRIs) inherent to quantum simulators introduce an extra layer of
complexity to the effective model realized, and it is vital to understand the
effect of LRIs on the efficacy of spatiotemporal quenches to implement them on
existing platforms.

For brevity, we focus our studies on spatiotemporal quenches in one specific
model—the long-range transverse field Ising (LR-TFI) model with ferromagnetic
power-law interactions and interaction amplitude $J(r)\sim 1/r^{\alpha}$, for
a range of $\alpha$. Previous work studying spatiotemporal quenches in the
short-range TFI and Heisenberg models [49, 1] whose low energy theories
exhibit Lorentz invariance, and thus a maximal speed of propagation of
information, and the presence of a linear lightcone that captures causality.
We investigate three experimentally motivated values of $\alpha$. We show that
for values of $\alpha=\\{3,6\\}$, where the critical dynamics are relativistic
[53], ‘Doppler-shift’ cooling persists and approximately the exact results for
free fermions if ultraviolet (UV) modes remain unexcited by the quench front.
Furthermore, our simulations reveal that there is a clear local minimum in the
energy density of excitations as a function of the velocity of the quench
front at $v\approx c$. These results highlight the efficacy of such
spatiotemporal quenches and underline the fact that optimal diabatic cooling
is obtained in the limit where the quench front velocity approaches the
emergent speed of light, $v\rightarrow c$.

Figure 1: Energy density during a spatiotemporal quench in the LR-TFI model
with $\alpha=6$ expressed in terms of the percentage of the critical spectral
bandwidth. The quench front acts as a source of excitations, populating modes
in a chiral way. At the end of the quench, the modes left in the wake of the
front are populated according to a red-shifted temperature leaving a “cold”
region of size $2\ell_{c}=Lc/v$ at energy $\epsilon_{c}$. Modes copropagating
with the front are populated according to a “hot” blue-shifted temperature
carrying energy $\epsilon_{h}$ and confined to a region of size
$2\ell_{h}=L(1-c/v)$.

For $\alpha=2$, the critical dynamics are non-relativistic; the protocol
accordingly loses its Doppler cooling effect. Instead, we focus on identifying
the relevant length scales that control correlations and excitations—we find
that both the usual QKZM length $\xi_{\text{KZ}}$ that governs the formation
of defect density in a homogeneous quench, and a length scale
$\xi_{\text{SP}}$ that governs relaxation of correlations near a domain wall
separating regions on either side of a phase transition, as introduced in Ref.
[49], are relevant. This generalized QKZM can be used to predict the behavior
of correlations and the energy density of excitations as the velocity of the
quench front is varied; we illustrate this using appropriate scaling
collapses. These generalized QKZM arguments in principle also apply to the
relativistic case but the scaling function itself has non-trivial behavior
that can only be predicted using an understanding of Doppler shifts.
Importantly, these Doppler shifts predict a local minimum in the excitation of
the system in the realitivistic case around $v\approx c$; such a minimum is
entirely absent for $\alpha=2$.

This manuscript is organized as follows. In Sec. II, we introduce the models
studied and their critical properties. In Sec. III, we show that the energy
density and correlation length at the end of the quench are qualitatively
consistent with Doppler cooling for $\alpha=3$ and $\alpha=6$. We compare our
findings to exact calculations computing quench dynamics in a system of free
relativistic fermions; these calculations are relegated to App. B for
readability of the main text. In Sec. IV, the case $\alpha=2$ is studied,
where using appropriate scaling collapses, we identify the length scales
governing excitations and correlation functions. We discuss the growth of
entanglement entropy during the quench at the end of each section. We conclude
with a summary of findings and potential future directions in Sec. V.

##  II Models studied and Quantum Kibble-Zurek scaling

We study spatiotemporal quenches in one-dimensional LR-TFI models with
Hamiltonians

$H=-J\left(\sum_{i<j}\frac{\sigma_{i}^{x}\sigma_{j}^{x}}{|i-j|^{\alpha}}+g_{c}\sum_{i}\sigma_{i}^{z}\right)-h\sum_{i}f_{i}(t)\sigma_{i}^{z},$
(1)

where $\sigma_{i}^{\mu}$ are the Pauli matrices, $J$ is the interaction
strength, $g_{c}$ is the critical transverse field and $h$ is the initial
gapping perturbation. In what follows, we set $J=1$ and let $h=4$ for all
three systems. The perturbation is quenched along smooth fronts moving at
velocity $v$ such that
$f_{i}(t)=\frac{1}{2}+\frac{1}{2}\tanh\\!\left[(\mathinner{\\!\left\lvert
x_{i}\right\rvert}-vt)/v\tau\right]$, where $\tau$ is the smoothing parameter.
The quench is started at time $t_{0}=-2\tau$ ensuring that
$f_{i}(t_{0})\approx 1$ at every site. At $t=\infty$, $f_{i}(t)=0\forall i$,
and the critical Hamiltonian is obtained. The quench time is halved by
starting the quench in the center of the chain. In this paper, we restrict the
study of quenches in chains of up to $N=256$ spins with all time-dependent
results presented for $N=128$.

The system is initialized in the ground state of the paramagnetic phase with
large $h$. The initial wavefunction is obtained using ITensor’s density matrix
renormalization group (DMRG) algorithm [54]. Concretely, the Hamiltonian is
represented as a matrix product operator (MPO) consisting of a sum of $k$
exponentially decaying Hamiltonians with different decay lengths that together
approximate power-law interactions [55, 56]. The error in the long-range
interaction amplitude at any site is restricted to $10^{-6}J$ by using
$k=3,6,8$ exponentials for $\alpha=6$, $3$ and $2$, respectively. The quality
of the approximation is verified by a scaling collapse of the energy gap for
all three systems—see App. A—where we find that the critical transverse fields
and critical exponents are consistent with previous numerical studies [57,
58]. Next, time evolution is carried out with the fourth-order time-dependent
variational principle [59, 60]. At every time step, the spin correlations, the
(von Neumann) entanglement entropy and the total energy with respect to the
critical Hamiltonian are calculated using standard matrix product state (MPS)
techniques [61].

We also compute the energy density; we define this locally over each bond
between sites $i,i+1$, as the expectation value of the operator

$\displaystyle h_{i}=$ $\displaystyle-\sum_{\text{odd }r\geq
1}\frac{1}{r^{\alpha}}\sigma^{x}_{i-\frac{r-1}{2}}\sigma^{x}_{i+\frac{r+1}{2}}$
$\displaystyle-\frac{1}{2}\sum_{\text{even }r\geq
2}\frac{1}{r^{\alpha}}\left(\sigma^{x}_{i-\frac{r}{2}}\sigma^{x}_{i+\frac{r}{2}}+\sigma^{x}_{i-\frac{r}{2}+1}\sigma^{x}_{i+\frac{r}{2}+1}\right)$
$\displaystyle-\frac{g_{c}}{2}\left(\sigma^{z}_{i}+\sigma^{z}_{i+1}\right)$
(2)

where $i\in[1,N-1]$. Here we used the critical value of the transverse field
as we are interested in finding the energy above the ground state of the
critical system. Note further that
$\sum_{i}h_{i}=H-\frac{g_{c}}{2}(\sigma^{z}_{1}+\sigma^{z}_{N}$; thus, the
local energy density defined in this way sums to the total energy besides
irrelevant boundary terms.

The choice of $\alpha$ studied in this work is motivated experimentally.
Dipole-dipole interactions and Van der Waals interactions between neutral
Rydberg atoms in optical traps naturally yield effective spin models with
$\alpha=3$ and $\alpha=6$ power-law interactions [62], respectively. Trapped
ion experiments can simulate effective spin models with phonon mediated long-
range interactions with $0<\alpha<3$. In practice, experiments are limited to
$\alpha\lesssim 1.5$ [9]. We note that in general one-dimensional spin models,
the tightest Lieb-Robinson bounds for $1<\alpha<2$ predict a logarithmic light
cone with a boundary $t\sim\log(r)$, polynomial light cones with $t\sim
r^{\kappa}$ for $2<\alpha\leq 3$ [63] and linear light cones $t\sim r$ for
$\alpha>3$ [64], with $\alpha=2$ and $\alpha=3$ being limiting cases of these
three regimes.

The LR-TFI model has a second order phase transition separating a
ferromagnetic phase for $g<g_{c}$ and paramagnetic phase for $g>g_{c}$. The
critical points of the Hamiltonians in Eq. (1) are found by performing a
scaling collapse of the energy gap calculated with DMRG as a function of the
transverse field, for various system sizes up to $N=192$; details are
presented in App. A. We find that the calculated critical fields $g_{c}$ and
critical exponents $\nu$ agree well with previous quantum Monte Carlo and DMRG
investigations [57, 58]. Moreover, renormalization group (RG) calculations of
the LR-TFI model with $\alpha=1+\sigma$ predict relativistic dynamics ($z=1$)
for $\sigma>7/4$. The critical theory for $\alpha\geq 3$ is understood to be
the same as the short range model, thus described by a free fermion theory.
For $2/3<\sigma<7/4$, the dynamical critical exponent is calculated in an
epsilon series expansion with $\epsilon=3\sigma/2-1$, giving
$z=\sigma/2+\rho(\sigma)\epsilon^{2}+O(\epsilon^{3})$ with
$\rho(\sigma)\approx 1/[24(1+\sigma^{2})]$. For $\alpha=2$, $z\approx 0.505$
and the critical dynamics are non-relativistic.

Near criticality, the equilibrium properties of the system become universal
and are described by a correlation length $\xi\sim|g-g_{c}|^{-\nu}$ where
$\nu$ is the critical exponent $\nu$ and $g-g_{c}$ measures the distance to
criticality. The corresponding energy gap is $\Delta\sim|g-g_{c}|^{z\nu}$,
where $z$ is the dynamical critical exponent. In homogeneous phase
transitions, where the parameter $g(t)$ is globally tuned to or across the
critical value $g_{c}$ at a rate $1/\tau$, the Kibble-Zurek mechanism
describes the typical length- and time-scales at which adiabatic evolution
breaks down [65, 66, 67]. Adiabatic evolution stops when the instantaneous
correlation length diverges at a rate faster than a threshold velocity. For a
general dispersion relation $\omega\sim q^{z}$, the group velocity of modes at
length scales $\xi$ (or $q\sim 1/\xi$) is given by $v_{q}\sim\xi^{1-z}$.
Solving $d\xi/dt=v_{q}$ gives $\xi_{\text{KZ}}\sim\tau^{\nu/(1+z\nu)}$ called
the healing length of the system. It describes the typical size of symmetry
broken regions caused by the excitations produced during the quench. The
threshold velocity can be estimated [49] as

$v^{*}\sim
v_{q}(q=1/\xi_{\text{KZ}})\sim\xi_{\text{KZ}}^{1-z}\sim\tau^{(1-z)\nu/(1+z\nu)}.$
(3)

In systems with $z=1$, the threshold velocity is a constant and given by the
maximum group velocity of quasiparticle excitations. The system evolves
adiabatically when the correlation length changes at a rate much slower than
this threshold velocity (the exact dependence of excitation energy on the
quench front velocity may exhibit important non-monotonicity as we discuss in
Sec. III).

In inhomogeneous phase transitions, the QKZM must be generalized to account
for the moving quench front [49]. At fixed time $T$, the quench front is at
location $x_{c}=vT$ and $g(x)\sim g_{c}+\frac{x-x_{c}}{v\tau}$ near the front.
Even at equilibrium, the correlation length diverges as $\xi(x)\sim
g(x)^{-\nu}$, which influences the behaviour of correlations near the front.
In a static system with an inhomogeneous perturbation $g(x)$, a length scale
$\xi_{\text{SP}}$ describing the decay of the order parameter to equilibrium
in the symmetry broken phase [68, 69] can be found by comparing the
instantaneous correlation length $\xi(x)$ with the distance $x-x_{c}$ to the
critical point on the quench front. This yields

$\xi_{\text{SP}}\sim(v\tau)^{\nu/(1+\nu)}.$ (4)

The correlation functions and the energy density in our system generally
depend on several length scales—the prominent ones being $\xi_{\text{KZ}}$ and
$\xi_{\text{SP}}$—and the spatial region in which the correlator is being
calculated. We attempt to identify these length scales in the various cases
studied and provide scaling functions where appropriate in the following
sections.

Figure 2:  (a) The total energy density for $\alpha=6$ at time $t_{q}$ is
shown. (b) The average energy density in the center of the spin chain is
shown. The averaging is done over the cold region $-R<x<R$, where
$R=\min(N/4,ct_{q})$ centered about $x=0$ and over a time interval $T=2$. The
red-dotted lines show the free fermion prediction with $c=2.036$ and the black
dotted line shows the energy density at the beginning of the quench.

##  III Doppler Cooling of Long-Range Models with $z\approx 1$

The heatwave picture developed in Ref. [50] summarizes the relativistic (or
Lorentz) cooling mechanism for superluminal quench fronts presented in this
section. In the heatwave picture, the population of modes excited by the
quench front is approximated by spatially segregated thermal distributions as
shown in Fig. 1. Modes copropagating with the front are confined to the region
$ct<\mathinner{\\!\left\lvert x\right\rvert}<vt$ and populated at a blue-
shifted temperature, while the counterpropagating modes occupying the region
$0<\mathinner{\\!\left\lvert x\right\rvert}<ct$ are populated at a red-shifted
temperature. For free relativistic fermions, as pertains to the short-range
limit of the critical TFI model, one can perform detailed calculation of these
mode populations; see App. B for details. The Pauli exclusion principle here
prevents the $\sim 1/k$ population of modes calculated in the case of
spatiotemporal quenches in free bosons [50, 1]. Instead, the population of
modes is given by

$N^{F}_{\theta}(k)\approx\frac{1}{2}\;\;\text{for}\;\;ck\ll\frac{m}{\gamma\eta(\theta)},$
(5)

where $\gamma=1/\sqrt{1-\beta^{2}}$ is the Lorentz factor with $\beta=c/v$,
$\eta(0)\equiv\eta=\sqrt{\frac{1+\beta}{1-\beta}}$ is the relativistic Doppler
factor associated with copropagating modes and $\eta(\pi)=1/\eta$ is
associated with counterpropagating modes. In the hot region, modes are
occupied up to $k\sim\frac{m}{c}\frac{\eta}{\gamma}$ while those in the cold
region are occupied up to $k\sim\frac{m}{c}\frac{1}{\eta\gamma}$. As
$v\rightarrow c^{+}$, the Doppler factor diverges and the cold region is left
completely unexcited. Integrating over momenta yields the result that the
energy density carried by right and left moving modes is $\sim
1/\left[\gamma^{2}\eta(\theta)^{2}\right]$.

Figure 3:  (a) The spin correlations for $\alpha=6$ are shown for $\tau=0.4$
and various velocities once the quench front reaches $x\approx 50$. The cold
region is identified with circles while the hot region is identified with
triangles. Crosses seen for $v=2$ represent sites that are still in the gapped
region locally (and thus have a decay length different from the cold and hot
regions altogether). (b) A scaling collapse of the spin correlations with
$\xi=\eta\gamma$ is shown for different $\tau$. The critical exponent of the
correlation $2\Delta=0.25$ is found by doing a scaling collapse of the
critical correlations (see appendix A). Note that for $\tau=\\{0.7,1\\}$, the
velocities $v<4.0$ are excluded from the scaling collapse as the correlation
length approaches $N/2$ and we are limited by finite size effects.

We note further that even in the nearest neighbor limit of the TFI model,
where an exact free fermion description holds, the mode dispersion deviates
from the relativistic form in the limit $k\rightarrow\pi$. These modes will
are not cooled with the above Doppler factors as their group velocity deviates
significantly from $c$. Thus, it is necessary to consider a smoothing
parameter $\tau$ (which is finite) that prevents UV modes from getting
excited—these are expected to be exponentially suppressed for
$ck\gg\tau\sim\mathcal{O}(1)$ which prevents heating at energies where the
dispersion relation deviates from the linear relation $\omega=ck$.

We begin by presenting results of spatiotemporal quenches for the LR-TFI
models with $\alpha=6$. A scaling collapse of the energy gap reveals that the
critical exponents exactly match those of free fermions ($z=1$, $\nu=1$) — see
App. A. The system is quenched to criticality with $g_{c}(N=128)=1.01$. Note
that $g_{c}(\infty)=1.031$ reported in App. A corresponds to the critical
value in the thermodynamic limit. The excitation energy density at the end of
the quench at time $t_{q}=N/2v+2\tau$ closely follows the theoretical
prediction for free fermions with speed of light $c=2.036$ as shown in Fig. 2.
For free fermions, (which should closely describe the critical properties of
the LR-TFI models with $\alpha\gtrsim 3$) we can compute exactly the energy
and spatial distribution of excitations due to the spatiotemporal quench. In
particular, the energy density $\epsilon_{c}$ in the cold region
$\mathinner{\\!\left\lvert x\right\rvert}<ct_{q}$ is found to be
$\epsilon_{c}\sim 1/\eta^{2}\gamma^{2}$ while that in the hot region
$ct_{q}<\mathinner{\\!\left\lvert x\right\rvert}<vt_{q}$ is
$\epsilon_{h}\sim\eta^{2}/\gamma^{2}$ — see App. B. The average energy density
over the entire length of the system must then be
$\epsilon_{\text{avg}}\approx\left(1-\frac{c}{v}\right)\epsilon_{h}+\frac{c}{v}\epsilon_{c}\sim\frac{1}{\gamma^{2}}$,
as shown in Fig. 2(a). Near $v=c^{+}$, the quench excites higher frequency
modes that are less efficiently cooled because of their nonlinear dispersion.
Increasing $\tau$ restores the cooling effect by adiabatically suppressing the
excitation of modes $ck\gg\tau^{-1}$ in the laboratory frame. This effect is
particularly strong in the cold region, where the energy density follows $\sim
1/\eta^{2}\gamma^{2}$ more closely as $\tau$ is increased as shown in Fig.
2(b).

Figure 4:  The growth of the von Neumann entanglement entropy during the
quench for $\alpha=6$ is shown for $\tau=1$. The red-dotted line shows linear
increase and the black-dotted line shows a logarithmic $\frac{1}{6}\log(ct)$
growth, as would be expected for the short range TFIM in its ground state over
a region of size $L=ct$. Inset: the entanglement entropy at the end of the
quench is shown for different $\tau$ as a function of front velocity. The
black-dotted line shows the entanglement calculated at criticality using DMRG.
Figure 5: The energy density during a spatiotemporal quench shown as the ratio
of $e(x,t)/e(t=0)$. (a) Quenches in systems with $\alpha=3$ show a persisting
heatwave picture. The white dotted line shows the light cone with $c$
estimated via a collapse of the spin correlations. The initial energy is
$e(t=0)=4.9\%$ of spectral bandwidth. (b) For $\alpha=2$, large front
velocities are necessary to see a separation of higher and lower energy
density regions. This is due to the existence of the dynamical threshold
velocity $v^{*}$ and is not a relativistic effect. Initially, $e(t=0)=1.8\%$
of spectral bandwidth. (c) Comparing $e(x,t_{q})/e(t=0)$ for different
$\alpha$. Figure 6:  (a) The average energy density for $\alpha=3$ at time
$t_{q}$ is shown. The averaging is done over the cold region $-R<x<R$, where
$R=\min(N/4,ct_{q})$ centered about $x=0$ and over a time interval $T=2$. The
red-dotted line shows the free fermion prediction for $c=3.705$ (as found in
appendix A and the black-dotted line shows that initial energy density. (b)
The spin correlations are collapsed using $\xi=\eta\gamma$. The best collapse
for $c=3.705$ is obtained for $2\Delta=0.17$, a value smaller than the
predicted $2\Delta=0.25$ (see appendix A). (c) The growth of von Neumann
entanglement entropy is shown. The red-dotted line shows linear increase while
the black-dotted shows a $1/6\log(t)$ increase.

The apparent minimum that can be seen near $v=c$ in Fig. 2(b) is due to
Doppler cooling in both limits of $v\rightarrow c$. For subluminal quenches,
the modes excited by the moving front continue interacting with it during the
entirety of the quench. The waves simply bounce back and forth between the
right and left propagating quench fronts. A classical solution to a moving
front is presented in App. C. It shows that in the limit of $v\rightarrow
c^{-}$, modes reflecting off of the moving front are red-shifted by
$1/\eta^{2}$ while no modes are transmitted. This leads to a $\sim 1/\eta^{2}$
cooling. Note that the Doppler cooling factor $\eta$ as defined in the
subluminal case is the same as in the superluminal case but with
$\beta\rightarrow 1/\beta$. Both lead to near-perfect red-shifts in the limit
$v\rightarrow c$.)

The heatwave picture motivated in Ref. [50] is also supported by the spin
correlations. The correlations decay exponentially with two length scales that
can be associated with the cold and hot regions as shown in Fig. 3(a). Using
the observation that the excitation energy scales as $\epsilon\sim\xi^{d+z}$
near criticality, the correlation length $\xi_{c}=\eta\gamma$ can be
identified in the cold region (corresponding to an energy density
$\epsilon_{c}$).

The correlator decays away from $x=0$ on two different length scales
corresponding to the cold and hot regions, as seen in Fig. 3(a). Additionally,
the correlation length increases as the quench front velocity $v$ approaches
$c$ in the cold region, while it decreases in the same limit in the hot
region, which agrees again with the heatwave picture. Quantitatively, we
examine the autocorrelator in the cold region using the following ansatz

$C_{x}(r<ct)=\langle\sigma_{0}^{x}\sigma_{r}^{x}\rangle=\xi_{c}^{-2\Delta}F_{C}\left(\frac{r}{\xi_{c}}\right),$
(6)

which identifies the scaling of the correlation length $\xi_{c}$ with the
quench velocity according to the heatwave picture. We note that the
autocorrelation function shows good scaling collapse with $\xi_{c}$ over
multiple quench front velocities, as seen in Fig. 3(b).

In the spatiotemporal quench we study, the entanglement entropy also appears
to increase slower than the expected linear growth in homogeneous quenches
[70, 71] as shown in Fig. 4. Optimal quench protocols with $v\rightarrow c$
and large $\tau$ show nearly logarithmic growth of the entanglement and
$S_{\text{vN}}$ at $t_{q}$ approaches the ground state value calculated at
criticality with DMRG, as shown in the inset of Fig. 4.

The heatwave picture persists in spatiotemporal quenches in the LR-TFI model
with $\alpha=3$, as shown in Fig. 5(a). The system is again initialized with
$h=4$ and quenched to criticality with $g_{c}(N=128)=1.38$. For this system,
RG calculations predict that $z=1$ and that critical dynamics are
relativistic. However, a fit of the minimum energy gap as a function of the
system size gives $z\approx 0.9$. The critical field and critical exponents
obtained from a collapse of the gap energy using $z=0.9$ closely matches the
results reported in Refs. [57, 58]. Despite this deviation from $z=1$, the
relativistic cooling effect persists as can be seen from the energy density in
the center of the chain at $t=t_{q}$ as shown in Fig. 6(a) where a clear
minimum is realized as a function of the quench front velocity, and the
appearance of hot and cold regions, as seen in Fig. 5(c). We note that in
principle a $\tau-$dependent threshold velocity $v^{*}$ as estimated by QKZM
arguments [Eq. (3)] can replace the critical velocity $c$ in a putative
heatwave picture that does not rely on relativistic Doppler cooling. However,
the velocity $v^{*}$ as understood in these arguments merely distinguishes
quenches which proceed adiabatically (for $v<v^{*}$) vs. those that proceed
non-adiabatically (for $v>v^{*}$). It cannot explain the _minimum_ observed in
the energy density as a function of the quench front velocity. The latter
requires a Doppler cooling interpretation and suggests the applicability of
such relativistic physics even in the $\alpha=3$ case.

We note further that the results are not in as close agreement with the free
fermion findings as for $\alpha=6$. The energy density in the cold region
appears to follow the free fermion result for velocities $v\gtrsim 5$ with the
critical velocity $c=3.705$ (calculated by fitting a light cone for
correlations arising from a local quench; see App. A) but fails to show good
agreement at smaller quench front velocities. Unlike the $\alpha=6$ case, the
interface between the hot and cold regions does not appear to coincide with
$\mathinner{\\!\left\lvert x\right\rvert}=ct$, as seen in Fig. 5(a), with the
velocity $c$ as identified above. This is putatively due to the longer range
of interactions which also smooth out the separation between hot and cold
regions as $\alpha$ is lowered, as seen in Fig. 5(c), and make identification
of a clear boundary between these regions challenging. A scaling collapse with
the correlation length $\xi=\eta\gamma$ is provided in Fig. 6(b). The best
collapse using $c=3.705$ reveals that $2\Delta=0.17$, which deviates slighlty
from the ground state value of $1/4$ calculated from a scaling collapse of the
spin correlations (see App. A). The discrepancy can arise from two sources—we
are not precisely at the ground state after these quenches, and the
correlation length can exceed the finite system size studied for a range of
velocities.

These results show in summary that the heatwave picture extends even to
$\alpha=3$, and the physics of Doppler cooling can be seen even though we
obtain a dynamical critical exponent that deviates slightly from $z=1$. The
entanglement entropy produced during the quench grows slower than linearly as
shown in Fig. 6(c) and grows nearly logarithmically for $v\lesssim 3$.
Together, these data show that the Doppler cooling effect in spatiotemporal
quenches continues to offer a strong advantage over homogeneous quenches in
preparing critical states of the LR-TFI model for $\alpha\gtrsim 3$ where
$z\approx 1$.

Figure 7:  (a) The total energy density for $\alpha=2$ at time $t_{q}$ is
shown to collapse according to Eq. 8. Inset The energy density at time $t_{q}$
is shown, with the black-dotted line indicating the initial energy density.
(b) The spin correlations near the quench front are collapsed for $v\leq 8$,
which shows that $\xi_{\text{SP}}$ is the dominant length scale. The quench
front is located at $x_{c}\approx 40$ in all cases. Inset The uncollapsed
correlations are shown, including $v=10$ and $15$. (c) The growth of von
Neumann entanglement entropy is shown. The red-dotted line shows linear
increase and the black-dotted line shows the entanglement entropy calculated
with DMRG at criticality.

##  IV Kibble-Zurek mechanism for $\alpha=2$

Increasing the range further to $\alpha=2$ leads to different dynamics. The
system is quenched to criticality with $g_{c}(N=128)=2.3285$. Higher and lower
energy regions, as shown in Fig. 5(b), are still visible in the quenches
although higher velocities are needed to observe this distinction. The
presence of hot and cold regions may be explained by the existence of the
dynamical threshold velocity $v^{*}$ in Eq. (3). As was the case with
$\alpha=3$, there is no sharp interface between these regions as shown in Fig.
5(c).

To analyze the properties of the system during the quench, we use the
dynamical critical exponent $z\approx 0.505$ calculated using RG in Ref. [53].
The critical theory is suggested to have dispersion modes with a dispersion
$\omega\sim q^{1/2}$, and Doppler-shift cooling arguments cannot be applied
here. Instead, we show that the dynamics are dominated by a competition
between the homogeneous-quench QKZM healing length
$\xi_{\text{KZ}}\sim\tau^{\nu/(1+z\nu)}$, and the static inhomogeneous length
scale $\xi_{\text{SP}}\sim(v\tau)^{\nu/(1+\nu)}$. The ratio of
$\xi_{\text{KZ}}$ and $\xi_{\text{SP}}$ in fact relates to the ratio of the
quench front velocity to the $\tau-$dependent threshold velocity of the
quench,

$\frac{\xi_{\text{SP}}}{\xi_{\text{KZ}}}\sim\left(\frac{v}{v^{*}}\right)^{\nu/(1+\nu)},$
(7)

where $v^{*}=\tau^{(1-z)\nu/(1+z\nu)}$ is the dynamical threshold velocity
identified in Eq. (3).

At fixed time $t^{*}$, a scaling ansatz for the spin correlations and
excitation energy density involving these two length scales yields

$\displaystyle\epsilon(t^{*})=\xi_{\text{SP}}^{-(d+z)}F_{e}\left(\frac{v}{v^{*}}\right),$
(8)
$\displaystyle\frac{C_{x}(x,t^{*})}{C_{x}(x_{c},t^{*})}=F_{C}\left(\frac{x-x_{c}}{\xi_{\text{SP}}},\frac{v}{v^{*}}\right),$
(9)

where $F_{e}$ and $F_{C}$ are unknown functions, and the ratio of the length
scales $\xi_{\text{KZ}}$ and $\xi_{\text{SP}}$ is considered inside the
scaling functions via the ratio $v/v^{*}$. Note that the scaling ansatz of
Eqs. (8,9) also works for $\alpha=\\{3,6\\}$ but only an exact calculation can
reveal the Doppler cooling effect that leads to non-monotonic cooling as a
function of $v$.

We begin by performing a scaling collapse of the total energy density at the
end of the quench according to Eq. (8) as shown in Fig. 7(a). The collapse
confirms that the dynamics are influenced by a non unique threshold velocity
$v^{*}$ determined by the quench rate $\tau$. We also note that there does not
appear to be any minimum in the energy density as a function of the quench
front velocity. The energy density simply decreases as this velocity is
lowered, in a marked difference from the result for $\alpha=3,6$. This further
confirms that the $\alpha=2$ system is not characterized by Doppler cooling. A
collapse of the spin correlations is done when the quench front reaches $x=40$
(at time $t=40/v$) according to Eq. (9), as shown in Fig. 7(b). The collapse
works reasonably well for $-10\xi_{\text{SP}}<x-x_{c}<5\xi_{\text{SP}}$ for
velocities $v<8$ and it shows that $\xi_{\text{SP}}$ is indeed the correct
length scale determining correlations close to the quench front. This length
scale also appears to play a role in determining the energy density over the
entire spin chain. For $v>8$, $\xi_{\text{SP}}$ is large enough at $x=40$ that
the collapse of correlations does not work well due to strong interaction with
the system boundaries, as can be seen in the inset of Fig. 7(b).

The von Neumann entanglement entropy also grows differently for $\alpha=2$
than for $\alpha=3,6$. The growth is approximately linear in time until the
quench front reaches the edges of the spin chain, where it saturates close to
the expected value of the entanglement entropy in the ground state of the
system (as found from DMRG), as can be seen in Fig. 7(c). However, the
entanglement entropy appears to show larger oscillations about this mean value
for larger velocities of the quench front, in agreement with the reduction in
energy density of excitations as the quench front velocity is lowered. [We
note in general that the entanglement entropy can be tricky to interpret as
low entanglement (as in a product state) and high entanglement (as found for
excited eigenstates) both imply heating.] The linear increase of
$S_{\text{vN}}$ is putatively due to the nearly logarithmic light cone
$t\sim\log r$ (see Ref. [63]) quickly spreading information over the entire
spin chain. It is an open question whether a faster protocol with a time-
dependent quench front velocity could be used to create low-energy states in
this system.

The above results largely imply that for $\alpha=2$, the front velocity serves
as another adiabatic parameter with optimal cooling achieved in the limit
$v\rightarrow 0$; a threshold velocity $v^{*}$ here can be used as a guideline
approximately demarcating adiabatic and non-adiabatic quenches.

##  V Conclusion

In this paper, we show that smooth spatiotemporal quenches can efficiently
prepare critical ground states of one dimensional LR-TFI models. We confirm
that for models with $z\approx 1$, when interactions $J(r)\sim 1/r^{\alpha}$
decay faster than the case $\alpha=3$, a smooth quench front moving along
$x=vt$ leaves a large section of the spin chain unexcited when $v\rightarrow
c$, the velocity of excitations in the critical system. This is evidenced by
the close resemblance between the energy density computed numerically for
these LR-TFI models and that calculated exactly for free relativistic
fermions. In particular, the excitation energy shows a local but pronounced
minimum as a function of the quench front velocity $v$ for $v\approx c$. For
general $v>c$, a heatwave picture emerges where one obtains spatially
separated hot and cold regions in the system populated by excitations
emanating from the quench front and either co- or counter-propagating with
respect to the moving quench front. Going from $\alpha=6$ to $\alpha=3$
smooths the separation between these hot and cool regions without
qualitatively impacting the results. Simple scaling relations accompanied by a
collapse of the spin correlations show that the correlations decay on the
diverging length scale $\xi_{c}\sim\eta\gamma$. The optimal quench protocol
also shows a nearly logarithmic growth of the von Neumann entanglement
entropy.

When the dynamical exponent $z$ deviates sufficiently from unity as for
$\alpha=2$ [$J(r)\sim 1/r^{2}$], one still obtains an approximate heatwave
picture of excitations in the system for large velocities of the quench front,
with areas of high and low energy density; a threshold velocity $v^{*}$
determined using QKZM arguments potentially serves the role played by the
critical velocity in the $z\approx 1$ case in separating hot and cold regions.
However, there is no local minimum in the energy density of excitations at a
function of the front velocity, which suggests an important departure from the
Doppler cooling picture prevalent in the above cases. Instead, we considered a
general scaling picture in this case—we find that energy density and spin
correlations are determined by the interplay between two length scales, one
corresponding to the size of broken symmetry clusters in a homogeneous quench
according to usual QKZM expectations, and another corresponding to the healing
length of correlations in a system with a spatially inhomogeneous gap with
$g<g_{c}$ on one side of the system and $g>g_{c}$ on the other side. The ratio
of these lengths is in fact related to the ratio of the velocity of the quench
front to a threshold velocity $v^{*}$ that controls adiabaticity of the
quench. The state prepared grows monotinically closer to the target critical
ground state as the velocity of the quench front is reduced.

This paper motivates the implementation of spatiotemporal quenches in one-
dimensional spin chains on modern quantum simulators. It remains to be shown
that Doppler cooling persists for two-dimensional systems with non-linear
dispersions. Spatiotemporal quenches are not expected to provide a cooling
advantage in all long-range systems. For power-law interactions with
$\alpha\leq 1$, Lieb-Robinson bounds do not exist and as such, there is no
light cone [63]. In that case, we could expect that the optimal protocol is
homogeneous as argued in Ref. [45] concerning the preparation of the critical
state of the fully connected ($\alpha=0$) TFI model. It however remains to be
understood whether a quench with a time-dependent front velocity $v(t)$ (see,
for instance, Ref. [2]) could be used to optimize the speed and efficacy of
ground state preparation for models with $z<1$ where some kind of light cone
exists.

###### Acknowledgements.

The authors acknowledge useful discussions with several previous collaborators
on related work. SB acknowledges support from an FRQNT graduate scholarship.
KA acknowledges support from the NSERC Discovery Grant and an INTRIQ team
grant from the FRQNT.

Figure 8:  The energy gap between the ground state and first excited state is
shown for spin chains of different length $N$. The insets show the collapsed
energy obtained with the indicated critical field $g_{c}$, critical exponent
$\nu$ and critical dynamical exponent $z$. (a) $\alpha=6$. (b) $\alpha=3$. (c)
$\alpha=2$. Figure 9: A scaling collapse of the spin correlations calculated
with DMRG at criticality for $\alpha={6,3,2}$. The scaling exponents $2\Delta$
are taken from the RG calculations in Ref. [53]. For $\alpha\geq 3$,
$2\Delta=0.25$ and for $\alpha=2$, $2\Delta\approx 0.505$. Figure 10: The
speed of light $c$ is estimated by doing a local perturbation of the ground
state of the LR-TFI model and time evolving with TDVP. The difference
$S_{\text{vN}}(x,t)-S_{\text{vN}}(x,t=0)$ evolves into a light cone following
the perturbation. The black-dotted lines represent the first instance of the
difference being larger than a cut-off of $0.1$. (a) In $\alpha=6$, $c\approx
2.036$. (b) In $\alpha=3$, $c\approx 3.705$. (c) In $\alpha=2$, the light-cone
is expected to be logarithmic and no attempt is made to determine the speed of
light.

##  Appendix A Properties of LR-TFI Models

The critical properties of the LR-TFI models presented in this work are
determined by doing a collapse of the energy gap between the ground state and
first excited state of the models calculated for different transverse fields
$h$ using DMRG. It also serves as a consistency check for the approximation of
power-law interacting Hamiltonian by a sum of exponentially decaying
Hamiltonians. For $\alpha=6$, we use the dynamical critical exponent $z=1$
based on Ref. [53]. The best collapse is obtained for the critical field
$g_{c}\approx 1.031$ and critical exponent $\nu=1$ as shown in Fig. 8(b),
closely matching previous investigations [57, 58].

For $\alpha=3$, the minimum gap at different $N$ is not consistent with $z=1$
as expected from RG calculations [53]. A fit of $\min_{h}\Delta=N^{-z}$ is
done to determine $z\approx 0.901$. Although this result differs from Ref.
[53], it is still consistent with the expectation that $\alpha=3$ is at the
limit between linear and polynomial Lieb-Robinson light-cones [64]. The best
collapse is obtained for $g_{c}\approx 1.424$ and $\nu\approx 1.050$ as shown
in Fig. 8(b), also consistent with previous studies [57, 58].

For $\alpha=2$, the dynamical critical exponent is chosen to be the corrected
mean-field value of $z\approx\sigma/2+\rho(\sigma)\epsilon^{2}$ calculated in
Ref. [53]. Here, $\sigma=1$, $\epsilon=3\sigma/2-1$ and $\rho(\sigma)\approx
1/[24(1+\sigma^{2})]$, giving $z\approx 0.505$. The best collapse gives
$g_{c}\approx 2.522$ and $\nu\approx 1.285$ as shown in Fig. 8(b), consistent
with the previously cited quantum Monte Carlo and DMRG studies [57, 58]. Out
of the three data sets shown in this section, $\alpha=2$ has the worst
collapse. We believe this is attributed to the large finite size effects
coming from using open boundary conditions with longer range interactions.

The scaling collapse of the spin correlations presented in Sec. III require
knowledge of the scaling dimensions of the correlator at criticality. We
calculate these critical correlations systems sizes of up to $N=256$ spins, as
shown in Fig. 9. We use the system size dependent critical fields $g_{c}(N)$,
which can be estimated by finding $\min_{h}\Delta$ at $N=128$. The collapse is
done using

$C_{x}(r)=N^{-2\Delta}F_{C}\left(\frac{r}{N}\right),$ (10)

where we assume the only relevant length at criticality for these systems is
the system size $N$. The collapse works well for the scaling dimensions
predicted by the RG calculations in Ref. [53]. For $\sigma>7/4$, RG predicts
that the scaling dimensions of $C_{x}$ is the same as the free fermion theory.
Thus, we use $2\Delta=0.25$ for $\alpha\geq 3$. For $2/3<\sigma<7/4$, the
scaling dimension is calculated using an epsilon series expansion giving
$2\Delta=1-\sigma/2+\rho(\sigma)\epsilon^{2}+O(\epsilon^{3})$. In the case of
$\alpha=2$, it gives $2\Delta\approx 0.505$, which collapses the data very
well.

The last property that we are interested in is the speed of “light” $c$ in
these systems. The following is inspired by previous work on the light-cone
spread of correlations in long-range interacting systems [72]. To include the
same finite size effects found in the spatiotemporal quenches, we restrict the
calculation to chains of $N=128$ spins. The critical ground states of the LR-
TFI models with $\alpha={2,3,6}$ are perturbed with the Pauli matrix
$\sigma^{z}$ at site $N/4$ (which we refer to as $x=0$), as shown in Fig. 10.
The wavefunction is evolved with the time-dependent variation principle
following the perturbation and the block von Neumann entanglement entropy is
calculated at each MPS bond at every time step. For $\alpha\geq 3$, the
difference $S_{\text{vN}}(x,t)-S_{\text{vN}}(x,t=0)$ evolves into a linear
light cone from which we can estimate the speed of light, as shown in Fig.
10(a-b). The estimate is done by finding when $S_{\text{vN}}(x,t)$ grows by a
certain amount (which we call the cut-off) compared to the ground state. The
cut-off is chosen to be $0.1$, a value that closely follows the light cone
boundary for $\alpha=3,6$. For $\alpha=2$, the light cone is expected to
spread as $t\sim\exp(3\sqrt{\log{r}})$ [63], such that no attempt at
calculating a speed of light was made.

##  Appendix B Formalism for Free Fermions and Solution to Instantaneous
Superluminal Quench

###  B.1 Problem Statement

The problem we would like to solve is specified by the following action and
commutation relations:

$\displaystyle
S=\int\\!dt\int^{L/2}_{-L/2}\\!dx\;\left[\frac{i}{2}\left(\bar{\psi}\gamma^{\alpha}\partial_{\alpha}\psi-\partial_{\alpha}\bar{\psi}\gamma^{\alpha}\psi\right)-m\Theta(x-v_{s}t)\bar{\psi}\psi\right]$
$\displaystyle\\{\psi_{a}(x,t),\psi_{b}^{\dagger}(x^{\prime},t)\\}=i\delta(x-x^{\prime})\delta_{a,b}\;\;\;\text{and}\;\;\;[\psi_{a}(x,t),\psi_{b}(x^{\prime},t)]=0.$
(11)

We set the speed of “light”, $c=1$. We will work in the Weyl basis wherein
$\gamma^{0}=\gamma_{x}$ and $\gamma^{1}=-i\sigma_{y}$, $\sigma_{x,y}$ being
Pauli matrices (see section 4.2.4 in Ref. [73], and more generally, Ref.
[74]), $\bar{\psi}=\psi^{\dagger}\gamma^{0}$, where $\psi_{a}(x,t)$ is the
field of the Weyl fermion with spinor index $a=1,2$ and
$\partial_{\mu}=(\partial_{t},\partial_{x})$. Note that we will use
superscript for the spinor index while the subscript will be used to
distinguish different solutions of the Dirac equation. The quench occurs
locally along a front that propagates towards the right at a fixed, supersonic
speed $v_{s}>1$. We define the inverse subsonic velocity $u_{s}\equiv
v^{-1}_{s}<1$.

While the Hamiltonian before and after the quench satisfies the usual discrete
symmetries associated with free relativistic fermions, imposing a boundary
necessitates the breaking of some of these symmetries, see Ref. [75]. We work
with ‘natural’ states that forgo parity (P) and charge-conjugate (C)
symmetries but keep time-reversal (T) and the combined CPT symmetry. These
states also satisfy the condition that the current
$c\psi^{\dagger}\sigma_{z}\psi$ is zero at the edges of the system. This is
affected with the following set of boundary conditions:
$\psi_{1}(L/2)=\psi_{2}(L/2)$, $\psi_{1}(-L/2)=-\psi_{2}(-L/2)$. One can check
for these conditions that $T\psi(x,t)=-\sigma_{x}\psi^{*}(x,-t)$ and
$CPT\psi(x,t)=-\sigma_{z}\psi(-x,-t)$ satisfy the same boundary conditions
while $P\psi(x,t)=\sigma_{x}\psi(-x,t)$ and
$C\psi(x,t)=\sigma_{z}\psi^{*}(x,t)$ do not.

###  B.2 Method of solution

####  B.2.1 General principle.

We work in the Heisenberg picture and describe the field operator prior to the
quench $(t<x/v_{s})$ by a mode-expansion in terms of the complete set of
solutions of the massive Dirac equation,
$i\gamma^{\mu}(\partial_{\mu}-m)\psi=0$. These are positive-frequency particle
(or ‘electron’) modes $v_{n}(x,t)$, and negative-frequency anti-particle (or
‘hole’) modes $\tilde{v}_{n}=Cv_{n}=\sigma_{z}v^{*}_{n}$, such that

$\psi(x,t<x/v_{s})=\sum_{n}\left[f_{n}v_{n}(x,t)+\tilde{f}^{\dagger}_{n}\tilde{v}_{n}(x,t)\right].$
(12)

The coefficients $f_{n}$ and $\tilde{f}_{n}$ satisfy the usual fermionic anti-
commutation relations: all operators anti-commute besides
$\\{f_{n},f^{\dagger}_{m}\\}=\delta_{n,m}$ and
$\\{\tilde{f}_{n},\tilde{f}^{\dagger}_{m}\\}=\delta_{n,m}$. The initial state
is defined via the relation $f_{n}\ket{0}=0$ and $\tilde{f}_{n}\ket{0}=0$ for
all $n$. Note that this amounts to setting the initial state as being the
vacuum of hole-like and particle-like excitations, which is the relevant case
for a critical system. As before, the above expansion is valid for all times
$t<x/v_{s}$ since this quench occurs on a space-like hypersurface.

After the quench, the field operator evolves according to the massless KG
equation and the mode expansion above is not valid for $t>x/v_{s}$. To find
correlations for subsequent times, we expand the massive modes in terms of the
massless modes. We reserve the notation $u_{n}(x,t)$ and $\tilde{u}_{n}(x,t)$
for the massless modes, and define

$\displaystyle v_{n}\big{|}_{x=v_{s}t}$
$\displaystyle=\sum_{m}\left[\alpha^{*}_{n,m}u_{m}+\beta_{n,m}\tilde{u}_{m}\right]\big{|}_{x=v_{s}t},$
$\displaystyle\tilde{v}_{n}\big{|}_{x=v_{s}t}$
$\displaystyle=\sum_{m}\left[\alpha_{n,m}\tilde{u}_{m}+\beta^{*}_{n,m}u_{m}\right]\big{|}_{x=v_{s}t},$
(13)

where $\alpha_{n,m}$ and $\beta_{n,m}$ are the Bogoliubov coefficients, and
the second equation follows from the first upon application of the charge-
conjugate operation. Then, the evolution of the field operator for times
$t>x/v_{s}$ can be described by the expansion

$\phi(x,t>x/v_{s})=\sum_{n}\left[\gamma_{n}(x,t)f_{n}+\tilde{\gamma}_{n}(x,t)\tilde{f}^{\dagger}_{n}\right],$
(14)

where

$\displaystyle\gamma_{n}(x,t)$
$\displaystyle=\sum_{m}\left[\alpha^{*}_{n,m}u_{m}(x,t)+\beta_{n,m}\tilde{u}_{m}(x,t)\right],$
$\displaystyle\tilde{\gamma}_{n}(x,t)$
$\displaystyle=\sum_{m}\left[\alpha_{n,m}\tilde{u}_{m}(x,t)+\beta^{*}_{n,m}u_{m}(x,t)\right].$
(15)

####  B.2.2 Dirac inner product and normalization of modes.

We use a coordinate-system invariant normalization scheme for the modes that
allows us to determine the Bogoliubov coefficients. We define the Dirac inner
product between two solutions $\psi_{a}$ and $\psi_{b}$ as

$\displaystyle(\psi_{a},\psi_{b})=\int dx\sqrt{g}n_{\mu}J^{\mu}_{(a,b)}(x),$
$\displaystyle\text{where}\quad
J^{\mu}_{(\psi_{a},\psi_{b})}=\bar{\psi}\gamma^{\mu}\phi.$ (16)

Here, the integral is over all space, $g$ is the determinant of the induced
metric on space-like coordinate, $\sqrt{g}dx$ is the covariant volume element,
$n^{\mu}$ is a future-directed time-like unit vector normal to the space-like
hypersurface and $J^{\mu}_{(a,b)}$ is the Dirac current. If $\psi_{a}$ and
$\psi_{b}$ satisfy the _same_ Dirac equation (massive or massless), then it is
easy to check that $\partial_{\mu}J^{\mu}_{(a,b)}=0$. Thus, the integral over
all space of the charge associated with the current $n_{\mu}J^{\mu}$ is
constant over time. Note that:

(a) If the modes $v_{n}$ and $\tilde{v}_{n}$ form a complete set of modes
according to the Dirac inner product, that is, $(v_{n},v_{m})=\delta_{n,m}$,
$(v_{n},\tilde{v}_{m})=0$, and the mode operators $f_{n}$ and $\tilde{f}_{n}$
satisfy the usual fermionic anti-commutation relations, then it can be shown
that the field operators (and its conjugate) satisfy the correct commutation
relations as described in Eq. (B.1).

(b) From its formulation in Eq. (B.2.2), it is explicit that the Dirac inner
product is invariant under transformation into a coordinate system which
admits a separation between time-like and space-like coordinates, that is, the
metric is of the form $ds^{2}=[N(x,t)]^{2}dt^{2}-g(x,t)dx^{2}$. Thus, the
normalization relations $(u_{n},u_{m})=\delta_{n,m}$,
$(v_{n},v_{m})=\delta_{n,m}$, etc. are invariant under such coordinate
transformations.

(c) The above two properties imply that under a Lorentz transformation of the
coordinates (without any change in the operators $f_{n}$, $\tilde{f}_{n}$),
the field operators continue to satisfy the commutation relations in Eq. (B.1)
in the _transformed_ coordinates, as appropriate for a relativistic field.

(d) The Dirac inner product has the symmetry that
$(\psi_{a},\psi_{b})=(\tilde{\psi}_{a},\tilde{\psi}_{b})$. Thus, both particle
and anti-particle modes follow the same normalization scheme
$(u_{n},u_{m})=(\tilde{u}_{n},\tilde{u}_{m})=\delta_{n,m}$.

####  B.2.3 Determination of $\alpha_{n,m}$ and $\beta_{n,m}$

To determine the coefficients $\alpha_{n,m}$ and $\beta_{n,m}$, we must
evaluate the Dirac inner product between modes along the curve $x=v_{s}t$. It
is useful to Lorentz-boost into a coordinate frame given by
$x^{\prime}=\gamma_{s}(x-u_{s}t)$, $t^{\prime}=\gamma_{s}(t-u_{s}x)$ with
$\gamma_{s}=1/\sqrt{1-u^{2}_{s}}$, as in this frame, the quench trajectory is
simply $t^{\prime}=0$.

The Dirac inner product evaluated at time $t^{\prime}=0$ in this frame reads

$(\psi_{a},\psi_{b})=\int_{-L/2\gamma_{s}}^{L/2\gamma_{s}}dx^{\prime}\;\psi_{a}^{\dagger}\psi_{b}\big{|}_{t^{\prime}=0}.$
(17)

Assuming Eqs. (B.2.1) hold at $t^{\prime}=0$, one may evaluate $(u_{n},v_{m})$
and $(u_{n},\tilde{v}_{m})$ to find

$\displaystyle(u_{n},v_{m})\big{|}_{t^{\prime}=0}$
$\displaystyle=\sum_{m}(u_{n},\alpha^{*}_{n,m}u_{m}+\beta_{n,m}\tilde{u}_{m})\big{|}_{t^{\prime}=0}=\alpha^{*}_{n,m},$
$\displaystyle(u_{n},\tilde{v}_{m})\big{|}_{t^{\prime}=0}$
$\displaystyle=\sum_{m}(u_{n},\alpha_{n,m}\tilde{u}_{m}+\beta^{*}_{n,m}u_{m})\big{|}_{t^{\prime}=0}=\beta^{*}_{n,m},$
(18)

where we used $(u_{n},\tilde{u}_{m})\big{|}_{t^{\prime}=0}=0$ and
$(u_{n},u_{m})\big{|}_{t^{\prime}=0}=\delta_{n,m}$. The above suggests that if
Eqs. (B.2.1) are simultaneously satisfiable, then the coefficients
$\alpha_{n,m}$ and $\beta_{n,m}$ must be given by Eq. (B.2.3). To confirm that
these are indeed the correct solutions, we can substitute these solutions into
Eqs. (B.2.1). Using the commutation relations on the field operators (as in
Eq. (B.1)) at $t^{\prime}=0$ directly confirms the validity of the result.

By the methods above, we may also show the inverse expansion at
$t^{\prime}=0$:

$\displaystyle u_{m}\big{|}_{x=v_{s}t}$
$\displaystyle=\sum_{n}\left[\alpha_{n,m}v_{n}+\beta_{n,m}\tilde{v}_{n}\right]\big{|}_{x=v_{s}t},$
$\displaystyle\tilde{u}_{m}\big{|}_{x=v_{s}t}$
$\displaystyle=\sum_{n}\left[\alpha^{*}_{n,m}\tilde{v}_{n}+\beta^{*}_{n,m}v_{n}\right]\big{|}_{x=v_{s}t}.$
(19)

Using these, one can easily prove that these fermionic Bogoliubov coefficients
have the following useful property:

$\displaystyle(u_{n},u_{m})$
$\displaystyle=\sum_{a}\left[\alpha^{*}_{a,n}\alpha_{a,m}+\beta^{*}_{a,n}\beta_{a,m}\right]=\delta_{n,m}.$
(20)

###  B.3 Solution of problem

####  B.3.1 Normalized modes.

We now provide details of the solution of the problem defined in Eq. (B.1).
The massive particle modes are defined as

$\displaystyle v_{\pm,k}$ $\displaystyle=\frac{1}{\sqrt{2L}}\left(v_{k}\pm
iv_{-k}\right)\;\text{for}\;\;k>0,$ $\displaystyle\text{where}\quad v_{k}$
$\displaystyle=\begin{pmatrix}\cos\\!\left(\theta_{k}/2\right)\\\
\sin\\!\left(\theta_{k}/2\right)\end{pmatrix}e^{-ikx+i\Omega_{k}t},$
$\displaystyle\tilde{v}_{-k}$
$\displaystyle=\begin{pmatrix}\sin\\!\left(\theta_{k}/2\right)\\\
\cos\\!\left(\theta_{k}/2\right)\end{pmatrix}e^{ikx+i\Omega_{k}t},$ (21)

while the anti-particle modes are defined as

$\displaystyle\tilde{v}_{\pm,k}$
$\displaystyle=\frac{1}{\sqrt{2L}}\left(\tilde{v}_{k}\pm
i\tilde{v}_{-k}\right)\;\text{for}\;\;k>0,$
$\displaystyle\text{where}\quad\tilde{v}_{k}$
$\displaystyle=\begin{pmatrix}\cos\\!\left(\theta_{k}/2\right)\\\
-\sin\\!\left(\theta_{k}/2\right)\end{pmatrix}e^{ikx-i\Omega_{k}t},$
$\displaystyle v_{-k}$
$\displaystyle=\begin{pmatrix}\sin\\!\left(\theta_{k}/2\right)\\\
-\cos\\!\left(\theta_{k}/2\right)\end{pmatrix}e^{-ikx-i\Omega_{k}t}.$ (22)

In the above,
$\cos(\theta_{k}/2)=\sqrt{\frac{1}{2}+\frac{1}{2}\frac{k}{\Omega_{k}}}$,
$\sin(\theta_{k}/2)=\sqrt{\frac{1}{2}-\frac{1}{2}\frac{k}{\Omega_{k}}}$, and
$\Omega_{k}=\sqrt{m^{2}+k^{2}}$. The modes satisfy the CPT symmetry conserving
boundary conditions for $kL=n\pi+\pi/2$, with $n\in[0,2,4,...)$ for modes
$v_{+,k}$ and $n\in[1,3,...)$ for modes $v_{-,k}$.

An analogous set of massless modes $u_{\pm,k>0}$ and $\tilde{u}_{\pm,k>0}$ can
be found by setting the mass to zero in the corresponding formulæ for the
massive modes. This is enforced by the substitutions
$\Omega_{k}\rightarrow\omega_{k}=\mathinner{\\!\left\lvert k\right\rvert}$,
$\cos(\theta_{k}/2)\rightarrow 1$, $\sin(\theta_{k}/2)\rightarrow 0$.

It is also useful to note the form of these modes in the Lorentz-boosted
frame. The coordinates and momenta are boosted in the usual way, with
$k^{\prime}x-\Omega_{k^{\prime}}t\rightarrow
k^{\prime}_{R}x^{\prime}-\Omega_{k^{\prime}_{R}}t^{\prime}$ and
$-k^{\prime}x-\Omega_{k^{\prime}}t\rightarrow-k^{\prime}_{L}x^{\prime}-\Omega_{k^{\prime}_{L}}t^{\prime}$.
The spinor part is transformed by multiplication with the matrix
$\Lambda\equiv-
ie^{\frac{\omega}{2}\cdot\frac{i}{4}[\gamma^{0},\gamma^{1}]}=\begin{pmatrix}1/\sqrt{\eta}&0\\\
0&\sqrt{\eta}\end{pmatrix}$. Note that $\omega=\text{tanh}^{-1}(-u_{s})$ is
the rapidity associated with the Lorentz boost and
$\eta=\sqrt{(1+u_{s})/(1-u_{s})}$ is the usual relativistic Doppler factor.

####  B.3.2 Bogoliubov coefficients.

The Bogoliubov coefficients can be evaluated by expressing these modes in the
Lorentz boosted coordinates and evaluating the Dirac inner product at time
$t^{\prime}=0$ as per Eq. (17) and Eqs. (B.2.3). The coefficients read

$\displaystyle\beta^{\epsilon,\epsilon^{\prime}}_{k,k^{\prime}}=(u_{\epsilon,k},\tilde{v}_{\epsilon^{\prime},k^{\prime}})^{*}=\frac{1}{2L}$
$\displaystyle\bigg{[}\frac{\cos\\!\left(\theta_{k^{\prime}}/2\right)}{\eta}F\\!\left(k_{R}+k^{\prime}_{R}\right)-i\epsilon^{\prime}\frac{\sin\\!\left(\theta_{k^{\prime}}/2\right)}{\eta}F\\!\left(k_{R}-k^{\prime}_{L}\right)$
$\displaystyle-i\epsilon\;\eta\;\sin\\!\left(\theta_{k^{\prime}}/2\right)F\\!\left(k^{\prime}_{R}-k_{L}\right)-\epsilon\epsilon^{\prime}\;\eta\;\cos\\!\left(\theta_{k^{\prime}}/2\right)F\\!\left(-k_{L}-k^{\prime}_{L}\right)\bigg{]}$
$\displaystyle\alpha^{\epsilon,\epsilon^{\prime}}_{k,k^{\prime}}=(u_{\epsilon,k},v_{\epsilon^{\prime},k^{\prime}})^{*}=\frac{1}{2L}$
$\displaystyle\bigg{[}\frac{\cos\\!\left(\theta_{k^{\prime}}/2\right)}{\eta}F\\!\left(k_{R}-k^{\prime}_{R}\right)-i\epsilon^{\prime}\frac{\sin\\!\left(\theta_{k^{\prime}}/2\right)}{\eta}F\\!\left(k_{R}+k^{\prime}_{L}\right)$
$\displaystyle+i\epsilon\;\eta\;\sin\\!\left(\theta_{k^{\prime}}/2\right)F\\!\left(-k^{\prime}_{R}-k_{L}\right)+\epsilon\epsilon^{\prime}\;\eta\;\cos\\!\left(\theta_{k^{\prime}}/2\right)F\\!\left(-k_{L}+k^{\prime}_{L}\right)\bigg{]}$
(23)

where
$F(x)=\frac{L}{\gamma_{s}}\text{sinc}\\!\left(\frac{xL}{2\gamma_{s}}\right)$
and the Doppler shifted momenta are given by

$\displaystyle k_{R}$
$\displaystyle=\gamma_{s}\left(k-u_{s}\omega_{k}\right)\in\big{[}\frac{\pi}{2L\eta},\infty\big{)},$
$\displaystyle k_{L}$
$\displaystyle=\gamma_{s}\left(k+u_{s}\omega_{k}\right)\in\big{[}\frac{\pi\eta}{2L},\infty\big{)},$
$\displaystyle k^{\prime}_{R}$
$\displaystyle=\gamma_{s}\left(k-u_{s}\Omega_{k}\right)\in[k^{\prime}_{0,-},\infty),$
$\displaystyle k^{\prime}_{L}$
$\displaystyle=\gamma_{s}\left(k+u_{s}\Omega_{k}\right)\in[k^{\prime}_{0,+},\infty),$
(24)

with $k^{\prime}_{0,\pm}=\gamma_{s}\frac{\pi}{2}(1\pm\sqrt{1+4m^{2}/\pi^{2}}$,
which correspond to frequencies

$\displaystyle\omega_{R}$
$\displaystyle=k_{R}=\gamma_{s}\left(\omega_{k}-u_{s}k\right),$
$\displaystyle\omega_{L}$
$\displaystyle=k_{L}=\gamma_{s}\left(\omega_{k}+u_{s}k\right),$
$\displaystyle\Omega_{R}$
$\displaystyle=\sqrt{k^{2}_{R}+m^{2}}=\gamma_{s}\left(\Omega_{k}-u_{s}k\right),$
$\displaystyle\Omega_{L}$
$\displaystyle=\sqrt{k^{2}_{L}+m^{2}}=\gamma_{s}\left(\Omega_{k}+u_{s}k\right).$
(25)

####  B.3.3 Infinite-size limit, chiral excitations populations.

We now work in the infinite-size limit and analyze the creation of excitations
from the vacuum. We look at the creation of massless anti-particles from the
massive particle modes sitting in the vacuum (and vice versa) since this
conversion precisely amounts to the excitation of the system about the vacuum
of the massless modes. Noting that the function $F(x)\rightarrow
2\pi\delta(x)$ in the limit $L\rightarrow\infty$, we find

$\displaystyle v_{\epsilon^{\prime},k^{\prime}}(t^{\prime}\\!=\\!0)$
$\displaystyle=\cos\\!\left(\frac{\theta_{k^{\prime}}}{2}\right)\left[\tilde{u}_{+,k}+\tilde{u}_{-,k}\right]\bigg{|}_{k=k^{-1}_{R}(-k^{\prime}_{R})}$
$\displaystyle-i\sin\\!\left(\frac{\theta_{k^{\prime}}}{2}\right)\left[\tilde{u}_{+,k}-\tilde{u}_{-,k}\right]\bigg{|}_{k=k^{-1}_{L}(k^{\prime}_{R})}$
$\displaystyle-i\epsilon^{\prime}\sin\\!\left(\frac{\theta_{k^{\prime}}}{2}\right)\left[\tilde{u}_{+,k}+\tilde{u}_{-,k}\right]\bigg{|}_{k=k^{-1}_{R}(k^{\prime}_{L})}$
$\displaystyle+\text{particle content}\propto u_{\pm,k}.$ (26)

The above is a direct result of the integration over the 4 momentum-conserving
delta-functions of the Bogoliubov coefficients and one of these terms,
$\propto\delta(-k_{L}-k^{\prime}_{L})$, does not contribute. Two of these
terms are associated with the production of right-movers
($\tilde{u}_{R,k}=\left(\tilde{u}_{+,k}+\tilde{u}_{-,k}\right)/\sqrt{2}$), and
one term is associated with left-movers
($\tilde{u}_{L,k}=-i\left(\tilde{u}_{+,k}-\tilde{u}_{-,k}\right)/\sqrt{2}$).
We focus on the first and second terms since the third term can be shown to be
continuously related to the first term, but carries a momentum $k>k_{0}$ while
the first carries momentum $k<k_{0}$, where $k_{0}=k^{-1}_{R}u_{s}\gamma_{s}m$
is of the order of the mass $m$.

We now evaluate the energy of the system after the quench. Note that the
Hamiltonian $H\equiv\psi^{\dagger}\cdot h\cdot\psi$, with
$h\equiv-i\sigma_{z}\partial_{x}$. For
$\psi(x,t>x/v_{s})=\sum_{\epsilon^{\prime},k^{\prime}}\left[\gamma_{\epsilon^{\prime},k^{\prime}}(x,t)f_{\epsilon^{\prime},k^{\prime}}+\tilde{\gamma}_{\epsilon^{\prime},k^{\prime}}(x,t)\tilde{f}^{\dagger}_{\epsilon^{\prime},k^{\prime}}\right]$,
and the state being a vacuum of operators $f_{\epsilon^{\prime},k^{\prime}}$
and $\tilde{f}_{\epsilon^{\prime},k^{\prime}}$, we find

$\displaystyle\left<H\right>$
$\displaystyle=\sum\left[\left(\alpha^{\epsilon_{1},\epsilon^{\prime}}_{k_{1},k^{\prime}}\right)^{*}\\!\left(\alpha^{\epsilon_{2},\epsilon^{\prime}}_{k_{2},k^{\prime}}\right)\\!\left(\tilde{u}^{\dagger}_{\epsilon_{1},k_{1}}\\!\cdot
h\cdot\tilde{u}_{\epsilon_{2},k_{2}}\right)\right]+\sum\left[\left(\beta^{\epsilon_{1},\epsilon^{\prime}}_{k_{1},k^{\prime}}\right)\\!\left(\beta^{\epsilon_{2},\epsilon^{\prime}}_{k_{2},k^{\prime}}\right)^{*}\\!\left(u^{\dagger}_{\epsilon_{1},k_{1}}\\!\cdot
h\cdot u_{\epsilon_{2},k_{2}}\right)\right]+\sim(\tilde{u}^{\dagger}\\!\cdot
u\text{ and }u^{\dagger}\\!\cdot\tilde{u})$
$\displaystyle=\sum\tilde{u}^{\dagger}_{\epsilon,k}\\!\cdot
h\cdot\tilde{u}_{\epsilon,k}-\sum\left(\sum\beta^{\epsilon,\epsilon^{\prime}}_{k,k^{\prime}}\tilde{u}_{\epsilon,k}\right)^{\dagger}\\!\cdot
h\cdot\left(\sum\beta^{\epsilon,\epsilon^{\prime}}_{k,k^{\prime}}\tilde{u}_{\epsilon,k}\right)+\sum\left(\sum\beta^{\epsilon,\epsilon^{\prime}}_{k,k^{\prime}}u_{\epsilon,k}\right)^{\dagger}\\!\cdot
h\cdot\left(\sum\beta^{\epsilon,\epsilon^{\prime}}_{k,k^{\prime}}u_{\epsilon,k}\right)$
$\displaystyle\approx\sum_{k}[1-N^{F}_{L}(k)](-k)\tilde{u}^{\dagger}_{L,k}\tilde{u}_{L,k}+\sum_{k}[1-N^{F}_{R}(k)](-k)\tilde{u}^{\dagger}_{R,k}\tilde{u}_{R,k}+\sum_{k}N^{F}_{L}(k)\;k\;u^{\dagger}_{L,k}u_{L,k}+\sum_{k}N^{F}_{R}(k)\;k\;u^{\dagger}_{R,k}u_{R,k}.$

In the above, all indices within any brackets are assumed to be summed over.
After the first equation, we neglect the oscillating (in time) terms of the
form $\tilde{u}^{\dagger}\cdot u$ and $u^{\dagger}\cdot\tilde{u}$ that are
expected to average out due to the integral over momenta. We also used Eq.
(20) to eliminate the $\alpha$ coefficients in favor of the $\beta$
coefficients. The last equation follows by substituting the result of Eq.
(B.3.3) into the second equation while neglecting time-dependent terms of the
form
$\tilde{u}^{\dagger}_{\epsilon,k}\tilde{u}_{\epsilon^{\prime},k^{\prime}}$ and
$u^{\dagger}_{\epsilon,k}u_{\epsilon^{\prime},k^{\prime}}$ that come with
$k\neq k^{\prime}$. The chiral populations read

$\displaystyle N^{F}_{L}(k)$
$\displaystyle=\frac{\Omega_{k_{L}}-\omega_{k_{L}}}{2\Omega_{k_{L}}},$
$\displaystyle N^{F}_{R}(k)$
$\displaystyle=\frac{\Omega_{k_{R}}-\omega_{k_{R}}}{2\Omega_{k_{R}}},$ (28)

where $k_{R}=k/\eta$ and $k_{L}=\eta k$.

Figure 11: (a) Population of fermions before the quench. The particle states
are unoccupied while the hole/anti-particle states are fully occupied. (b) The
population of hole and particle states is changed after the quench for momenta
$k<m/\eta$ for left-moving modes and $k<\eta m$ for right-moving modes. In the
limit $v_{s}\rightarrow 1^{+}$, $\eta\rightarrow\infty$, and we see that the
left-moving modes are left entirely unexcited.

There is a cooling effect for fermions because the population starts decaying
at a Doppler-shifted energy scale: $N^{F}_{L}(k)\sim 1/2$ for $k\ll m/\eta$
while $N^{F}_{R}(k)\sim 1/2$ for $k\ll\eta m$, and both decrease as $1/k^{2}$
for larger $k$. The population of fermions before and after the quench is
illustrated in Fig. (11). The factor of $1/2$ occurs due to spinor overlap
between the massive and massless modes: at $k=0^{+}$, the massless modes have
a spinor wave-function $(1,0)^{T}$ or $(0,1)^{T}$, but the massive modes have
a wave-function $\frac{1}{\sqrt{2}}(1,\pm 1)^{T}$. Thus, the overlap cannot
exceed $1/2$ due to Pauli exclusion. Another important difference is that the
above result has a UV singularity $\sim 1/k^{2}$. Thus, the result corresponds
to a UV singularity even in $d=1$. This UV singularity can be eliminated by
adding a time-scale to the quench.

###  B.4 Quench with a finite time-scale

In this section, we discuss how the $1/k^{2}$ UV singularity is removed by
adding a time-scale to the quench. This problem has been analyzed in the
context of particle production in inflationary cosmology [76]. There, the
metric undergoes a scale change that is equivalent to a scaling of the mass as
it is the only term that breaks the conformal invariance of the theory. For
the sake of completeness, we note the method used there and quote the result
relevant for our purposes.

We again study the case where the quench occurs uniformly in all space and
impose periodic boundary conditions. A generic solution to the time-dependent
Dirac equation of motion is then given by

$\displaystyle\chi_{k}(t)=\frac{1}{\sqrt{L}}\big{[}$
$\displaystyle\alpha_{k}(t)\begin{pmatrix}\cos(\theta_{k}(t)/2)\\\
\sin(\theta_{k}(t)/2)\end{pmatrix}e^{-i\int^{t}\omega_{k}(t^{\prime})dt^{\prime}}$
$\displaystyle+$
$\displaystyle\beta_{k}(t)\begin{pmatrix}\sin(\theta_{k}(t)/2)\\\
\cos(\theta_{k}(t)/2)\end{pmatrix}e^{\int^{t}\omega_{k}(t^{\prime})dt^{\prime}}\big{]},$
(29)

where $\alpha_{k}(-\infty)=1$, $\beta_{k}(-\infty)=0$,
$\omega_{k}(t)=\sqrt{m^{2}(t)+\mathinner{\\!\left\lvert k\right\rvert}^{2}}$
and the angle $\theta_{k}(t)$ is decided by the instantaneous frequency of the
mode $\omega_{k}(t)$. Note that, $\chi_{k}(t\rightarrow-\infty)$ reduces to
the massive mode solution $v_{k}$, while in the $t\rightarrow\infty$ limit it
is a linear combination of massless particle and anti-particle solutions with
momentum $k$. We are interested in
$\mathinner{\\!\left\lvert\beta_{k}(t=+\infty)\right\rvert}^{2}$ which is the
population of excitations at momentum $k$ after the quench is over.

Plugging in the ansatz of Eq. (B.4) into the time-dependent Dirac equation of
motion, we find

$\displaystyle\frac{d\alpha_{k}}{dt}$
$\displaystyle=-\beta_{k}\frac{k\;dm/dt}{2\omega^{2}_{k}}e^{2i\int^{t}\omega_{k}(t^{\prime})dt^{\prime}},$
$\displaystyle\frac{d\beta_{k}}{dt}$
$\displaystyle=\alpha_{k}\frac{k\;dm/dt}{2\omega^{2}_{k}}e^{-2i\int^{t}\omega_{\bm{k}}(t^{\prime})dt^{\prime}}.$
(30)

We assume $\beta_{k}(t)\ll 1$ and $\alpha_{k}(t)\approx 1$ (justified a
posteriori) and solve for
$\mathinner{\\!\left\lvert\beta_{k}\right\rvert}^{2}$. We now use the result
from Ref. [76] for this integral: it is approximated using the steepest
descent method and is a reasonable approximation when
$\tau^{-1}\lesssim\omega_{k}(t)$ — thus it is valid for $\tau^{-1}\lesssim m$
and $k\gtrsim m$. The result for the mass
$m(t)=mf(-t/\tau)=\frac{1}{2}+\frac{1}{2}\text{tanh}(-t/\tau)$ is

$N_{k}\underset{k\gg\tau^{-1},m}{\approx}\mathinner{\\!\left\lvert\beta_{k}\right\rvert}^{2}\approx
e^{-2m\tau-2\frac{k^{2}}{m}\tau}.$ (31)

Thus, the excitation of modes with momentum $k\gtrsim\tau^{-1}$ is suppressed
exponentially. The case where the quench occurs non-uniformly, via a space-
and time-dependent mass $m(x,t)=mf[(x-v_{s})/(v_{s}\tau)]$ cannot be analyzed
exactly due to the lack of momentum conservation. While one can boost to a
frame in which the quench does occur uniformly, the boundaries in this frame
can provide momentum kicks, making an exact analysis difficult. Nevertheless,
we expect that the analysis above should remain valid for modes with momenta
$\gg 1/L$ that are not particularly sensitive to the boundaries.

In the boosted frame, the quench occurs as
$m(x^{\prime},t^{\prime})=mf\left[-\frac{t^{\prime}}{\tau\gamma_{s}}\right]$.
Thus, we expect the inverse of the time-scale in the boosted frame,
$\tau^{\prime-1}=\tau^{-1}/\gamma_{s}$ to become the relevant energy scale
above which excitations are suppressed _in the boosted frame_. In the
laboratory frame, this implies a population

$\displaystyle N^{F}_{\theta}(k)$
$\displaystyle\approx\frac{1}{2}\;\;\text{for}\;\;k\ll\frac{m}{\gamma_{s}\eta(\theta)}$
(32)

and decaying exponentially in the opposite limits. To remind the reader, for
left-moving modes, $\theta=0$, $\eta(\theta)=\eta$, and for right-moving
modes, $\theta=\pi$, $\eta(\theta)=1/\eta$. Thus, the energy density of the
left- and right-moving modes is given by

$\epsilon_{\theta}(\tau)\bigg{|}_{\tau^{-1}\approx
m}\propto\int^{\frac{m}{\gamma_{s}\eta(\theta)}}\frac{k}{2}dk\propto\frac{m}{L_{m}}\frac{1}{(\gamma_{s}\eta(\theta))^{2}},$
(33)

where $L_{m}=m/c$ has dimensions of length.

##  Appendix C Subluminal Cooling

In spatiotemporal quenches in the LR-TFI models with $z=1$, the energy density
shows a minimum when the speed $v$ of the front approaches the speed of light.
The superluminal case is studied in appendix B. Here, we offer a classical
argument for the existence of this minimum due to relativistic effects in the
subluminal case. In the following, we assume that $v<c$. We note here as well
that in the limit $v\rightarrow 0$, the energy density will again be
suppressed, this time because it corresponds to the adiabatic limit.

The modes excited by a subluminal moving quench front interact with both
fronts throughout the entirety of the quench—the excitations released from,
say, the right-moving quench front will eventually bounce off the left-moving
quench front and interact again with the right-moving front. Classically, the
incoming wave will be transmitted through and reflected by the front. The
moving front forces moving boundary conditions that change the transmitted and
reflected waves frequencies and momenta. The following is based on previous
work on front induced wavepacket engineering [77, 78]. The fields on both
sides of the moving front are given by

$\displaystyle\phi_{-}(x,t)$
$\displaystyle=a_{\text{in}}e^{-i(\omega_{\text{in}}t-k_{\text{in}}x)}+a_{r}e^{-i(\omega_{r}t+k_{r}x)}$
(34) $\displaystyle\phi_{+}(x,t)$
$\displaystyle=a_{t}e^{-i(\omega_{t}t-k_{t}x)}$ (35)

where $a_{i}$ are the amplitudes, $\omega_{i}$ the frequencies and $k_{i}$ the
wavevectors of the incoming, reflected, and transmitted waves. The boundary
conditions at $x=vt$ are $D\phi_{+}=D\phi_{-}$ where
$D=\\{1,\partial_{t},\partial_{x}\\}$ leading to the phase matching conditions

$\displaystyle\frac{\omega_{r}}{\omega_{\text{in}}}$
$\displaystyle=\frac{1-v/V_{\text{in}}}{1+v/V_{r}}$ (36)
$\displaystyle\frac{\omega_{t}}{\omega_{\text{in}}}$
$\displaystyle=\frac{1-v/V_{\text{in}}}{1-v/V_{t}}$ (37)

where $V_{i}=\omega_{i}/k_{i}$ is the phase velocity of the different waves.

The modes in the $x<vt$ region are massless and have dispersion $\omega=ck$.
The phase velocity is simply $V=c$ such that Eq.(36) becomes

$\frac{\omega_{r}}{\omega_{\text{in}}}=\frac{1-\beta}{1+\beta}=\frac{1}{\eta^{2}},$
(38)

where $\beta=v/c$, which we recognize as a Doppler red-shifting of the
frequency. Therefore, the energy of reflected modes is red-shifted leading to
a strong cooling effect as $v\rightarrow c^{-}$. We note that the Doppler
shift factor employed here is different (with $\beta$ inversed to be precise)
from the one employed in the superluminal case, but the net result is the
same—cooling is strongest close to $v\rightarrow c$.

The modes transmitted into massive region $x>vt$ acquire a mass such that
their dispersion becomes $\omega=c\sqrt{m^{2}+k^{2}}$. The phase velocity is
$V_{t}=c/\sqrt{1-m^{2}/\omega_{t}^{2}}$ and is now frequency dependent.
Plugging it into Eq.(37) and solving for $\omega_{t}$ gives

$\frac{\omega_{t}}{\omega_{\text{in}}}=\frac{1+n^{\prime}\beta}{1+\beta},$
(39)

where

$n^{\prime}=\sqrt{1-\frac{\eta^{2}m^{2}}{\omega_{\text{in}}^{2}}}$ (40)

is used to define a cut-off frequency. Indeed, only modes with $\eta
m<\omega_{\text{in}}$ can be transmitted into the massive region. In the limit
$v\rightarrow c^{-}$, $\eta$ diverges and no modes are transmitted. Therefore,
the energy density depends only on the modes excited and reflected by the
moving quench front. Having shown that the latter is red-shifted, the energy
density is dominated by the modes excited by the subluminal moving quench
front, a subject of future investigation.

## References

  * Agarwal  _et al._ [2018] K. Agarwal, R. N. Bhatt, and S. L. Sondhi, Fast Preparation of Critical Ground States Using Superluminal Fronts, Phys. Rev. Lett. 120, 210604 (2018).
  * Mitra  _et al._ [2019] P. Mitra, M. Ippoliti, R. N. Bhatt, S. L. Sondhi, and K. Agarwal, Cooling arbitrary near-critical systems using hyperbolic quenches, Physical Review B 99, 104308 (2019).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Reviews of Modern Physics 80, 885 (2008).
  * Bloch  _et al._ [2012] I. Bloch, J. Dalibard, and S. Nascimbène, Quantum simulations with ultracold quantum gases, Nature Physics 8, 267 (2012).
  * Gross and Bloch [2017] C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science 357, 995 (2017).
  * Browaeys and Lahaye [2020] A. Browaeys and T. Lahaye, Many-body physics with individually controlled Rydberg atoms, Nature Physics 16, 132 (2020).
  * Blatt and Roos [2012] R. Blatt and C. F. Roos, Quantum simulations with trapped ions, Nat. Phys. 8, 277 (2012).
  * Gärttner  _et al._ [2017] M. Gärttner, J. G. Bohnet, A. Safavi-Naini, M. L. Wall, J. J. Bollinger, and A. M. Rey, Measuring out-of-time-order correlations and multiple quantum spectra in a trapped-ion quantum magnet, Nature Physics 13, 781 (2017).
  * Monroe  _et al._ [2021] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Programmable quantum simulations of spin systems with trapped ions, Rev. Mod. Phys. 93, 025001 (2021).
  * Agarwal and Martin [2020] K. Agarwal and I. Martin, Dynamical enhancement of symmetries in many-body systems, Phys. Rev. Lett. 125, 080602 (2020).
  * Keesling  _et al._ [2019] A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, P. Zoller, M. Endres, M. Greiner, V. Vuletic, and M. D. Lukin, Quantum Kibble-Zurek mechanism and critical dynamics on a programmable Rydberg simulator, Nature 568, 207 (2019).
  * Jepsen  _et al._ [2021] P. N. Jepsen, W. W. Ho, J. Amato-Grill, I. Dimitrova, E. Demler, and W. Ketterle, Transverse Spin Dynamics in the Anisotropic Heisenberg Model Realized with Ultracold Atoms, Phys. Rev. X 11, 041054 (2021).
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Quantum simulation of 2D antiferromagnets with hundreds of Rydberg atoms, Nature 595, 233 (2021).
  * Ebadi  _et al._ [2021] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletic, and M. D. Lukin, Quantum phases of matter on a 256-atom programmable quantum simulator, Nature 595, 227 (2021).
  * Scholl  _et al._ [2022] P. Scholl, H. J. Williams, G. Bornet, F. Wallner, D. Barredo, L. Henriet, A. Signoles, C. Hainaut, T. Franz, S. Geier, A. Tebben, A. Salzinger, G. Zürn, T. Lahaye, M. Weidemüller, and A. Browaeys, Microwave Engineering of Programmable $XXZ$ Hamiltonians in Arrays of Rydberg Atoms, PRX Quantum 3, 020303 (2022).
  * Tarruell  _et al._ [2012] L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, and T. Esslinger, Creating, moving and merging Dirac points with a Fermi gas in a tunable honeycomb lattice, Nature 483, 302 (2012).
  * Grusdt  _et al._ [2013] F. Grusdt, M. Höning, and M. Fleischhauer, Topological Edge States in the One-Dimensional Superlattice Bose-Hubbard Model, Physical Review Letters 110, 260405 (2013).
  * de Leseleuc  _et al._ [2019] S. de Leseleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Buchler, T. Lahaye, and A. Browaeys, Observation of a symmetry-protected topological phase of interacting bosons with Rydberg atoms, Science 365, 775 (2019).
  * Wintersperger  _et al._ [2020] K. Wintersperger, C. Braun, F. N. Ünal, A. Eckardt, M. D. Liberto, N. Goldman, I. Bloch, and M. Aidelsburger, Realization of an anomalous Floquet topological system with ultracold atoms, Nature Physics 16, 1058 (2020).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021).
  * Pal and Huse [2010] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010).
  * Alet and Laflorencie [2018] F. Alet and N. Laflorencie, Many-body localization: An introduction and selected topics, Comptes Rendus Physique 19, 498 (2018).
  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Observation of many-body localization of interacting fermions in a quasirandom optical lattice, Science 349, 842 (2015).
  * Choi  _et al._ [2016] J. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Exploring the many-body localization transition in two dimensions, Science 352, 1547 (2016).
  * Smith  _et al._ [2016] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Many-body localization in a quantum simulator with programmable random disorder, Nature Physics 12, 907 (2016).
  * Morong  _et al._ [2021] W. Morong, F. Liu, P. Becker, K. S. Collins, L. Feng, A. Kyprianidis, G. Pagano, T. You, A. V. Gorshkov, and C. Monroe, Observation of Stark many-body localization without disorder, Nature 599, 393 (2021).
  * Zohar  _et al._ [2016] E. Zohar, J. I. Cirac, and B. Reznik, Quantum simulations of lattice gauge theories using ultracold atoms in optical lattices, Reports on Progress in Physics 79, 014401 (2016).
  * Hauke  _et al._ [2013] P. Hauke, D. Marcos, M. Dalmonte, and P. Zoller, Quantum Simulation of a Lattice Schwinger Model in a Chain of Trapped Ions, Physical Review X 3, 041018 (2013), pRX.
  * Martinez  _et al._ [2016] E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl, P. Hauke, M. Dalmonte, T. Monz, P. Zoller, and R. Blatt, Real-time dynamics of lattice gauge theories with a few-qubit quantum computer, Nature 534, 516 (2016).
  * Bañuls  _et al._ [2020] M. C. Bañuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, C. A. Muschik, B. Reznik, E. Rico, L. Tagliacozzo, K. Van Acoleyen, F. Verstraete, U.-J. Wiese, M. Wingate, J. Zakrzewski, and P. Zoller, Simulating lattice gauge theories within quantum technologies, The European Physical Journal D 74, 165 (2020).
  * Choi  _et al._ [2019] J. Choi, H. Zhou, S. Choi, R. Landig, W. W. Ho, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, D. A. Abanin, and M. D. Lukin, Probing Quantum Thermalization of a Disordered Dipolar Spin Ensemble with Discrete Time-Crystalline Order, Phys. Rev. Lett. 122, 043603 (2019).
  * Bluvstein  _et al._ [2021] D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho, S. Choi, M. Serbyn, M. Greiner, V. Vuletic, and M. D. Lukin, Controlling quantum many-body dynamics in driven Rydberg atom arrays, Science 371, 1355 (2021).
  * Viermann  _et al._ [2022] C. Viermann, M. Sparn, N. Liebster, M. Hans, E. Kath, A. Parra-López, M. Tolosa-Simeón, N. Sánchez-Kuntz, T. Haas, H. Strobel, S. Floerchinger, and M. K. Oberthaler, Quantum field simulator for dynamics in curved spacetime, Nature 611, 260 (2022).
  * Kasper  _et al._ [2021] V. Kasper, D. González-Cuadra, A. Hegde, A. Xia, A. Dauphin, F. Huber, E. Tiemann, M. Lewenstein, F. Jendrzejewski, and P. Hauke, Universal quantum computation and quantum error correction with ultracold atomic mixtures, Quantum Science and Technology 7, 015008 (2021).
  * Ebadi  _et al._ [2022] S. Ebadi, A. Keesling, M. Cain, T. T. Wang, H. Levine, D. Bluvstein, G. Semeghini, A. Omran, J.-G. Liu, R. Samajdar, X.-Z. Luo, B. Nash, X. Gao, B. Barak, E. Farhi, S. Sachdev, N. Gemelke, L. Zhou, S. Choi, H. Pichler, S.-T. Wang, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum optimization of maximum independent set using Rydberg atom arrays, Science 376, 1209 (2022).
  * Mazurenko  _et al._ [2017] A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, A cold-atom Fermi–Hubbard antiferromagnet, Nature 545, 462 (2017).
  * Linke  _et al._ [2018] N. M. Linke, S. Johri, C. Figgatt, K. A. Landsman, A. Y. Matsuura, and C. Monroe, Measuring the Rényi entropy of a two-site Fermi-Hubbard model on a trapped ion quantum computer, Phys. Rev. A 98, 052334 (2018).
  * Tarruell and Sanchez-Palencia [2018] L. Tarruell and L. Sanchez-Palencia, Quantum simulation of the Hubbard model with ultracold fermions in optical lattices, C. R. Phys. 19, 365 (2018).
  * Bohrdt  _et al._ [2021] A. Bohrdt, L. Homeier, C. Reinmoser, E. Demler, and F. Grusdt, Exploration of doped quantum magnets with ultracold atoms, Annals of Physics 435, 168651 (2021), special issue on Philip W. Anderson.
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Adiabatic quantum computation, Rev. Mod. Phys. 90, 015002 (2018).
  * del Campo  _et al._ [2012] A. del Campo, M. M. Rams, and W. H. Zurek, Assisted Finite-Rate Adiabatic Passage Across a Quantum Critical Point: Exact Solution for the Quantum Ising Model, Phys. Rev. Lett. 109, 115703 (2012).
  * Damski [2014] B. Damski, Counterdiabatic driving of the quantum Ising model, J. Stat. Mech. 2014, P12019 (2014).
  * Sels and Polkovnikov [2017] D. Sels and A. Polkovnikov, Minimizing irreversible losses in quantum systems by local counterdiabatic driving, Proc. Natl. Acad. Sci. 114, E3909 (2017).
  * Pichler  _et al._ [2018] H. Pichler, S.-T. Wang, L. Zhou, S. Choi, and M. D. Lukin (2018), arXiv:1808.10816 [quant-ph].
  * Ho  _et al._ [2019] W. W. Ho, C. Jonay, and T. H. Hsieh, Ultrafast variational simulation of nontrivial quantum states with long-range interactions, Phys. Rev. A 99, 052332 (2019).
  * Pagano  _et al._ [2020] G. Pagano, A. Bapat, P. Becker, K. S. Collins, A. De, P. W. Hess, H. B. Kaplan, A. Kyprianidis, W. L. Tan, C. Baldwin, L. T. Brady, A. Deshpande, F. Liu, S. Jordan, A. V. Gorshkov, and C. Monroe, Quantum approximate optimization of the long-range Ising model with a trapped-ion quantum simulator, Proc. Natl. Acad. Sci. USA 117, 25396 (2020).
  * Ho and Zhou [2009] T.-L. Ho and Q. Zhou (2009), arXiv:0911.5506 [cond-mat.quant-gas].
  * Zaletel  _et al._ [2021] M. P. Zaletel, A. Kaufman, D. M. Stamper-Kurn, and N. Y. Yao, Preparation of Low Entropy Correlated Many-Body States via Conformal Cooling Quenches, Phys. Rev. Lett. 126, 103401 (2021).
  * Dziarmaga and Rams [2010] J. Dziarmaga and M. M. Rams, Dynamics of an inhomogeneous quantum phase transition, New J. Phys. 12, 055007 (2010).
  * Agarwal  _et al._ [2017] K. Agarwal, E. G. Dalla Torre, J. Schmiedmayer, and E. Demler, Quantum heat waves in a one-dimensional condensate, Phys. Rev. B 95, 195157 (2017).
  * Auerbach [2012] A. Auerbach,  _Interacting Electrons and Quantum Magnetism_ (Springer New York, 2012).
  * Giamarchi [2003] T. Giamarchi,  _Quantum physics in one dimension_ , Vol. 121 (Clarendon press, 2003).
  * Maghrebi  _et al._ [2016] M. F. Maghrebi, Z.-X. Gong, M. Foss-Feig, and A. V. Gorshkov, Causality and quantum criticality in long-range lattice models, Phys. Rev. B 93, 125128 (2016).
  * Fishman  _et al._ [2022] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022).
  * Crosswhite  _et al._ [2008] G. M. Crosswhite, A. C. Doherty, and G. Vidal, Applying matrix product operators to model systems with long-range interactions, Phys. Rev. B 78, 035116 (2008).
  * Pirvu  _et al._ [2010] B. Pirvu, V. Murg, J. I. Cirac, and F. Verstraete, Matrix product operator representations, New J. Phys 12, 025012 (2010).
  * Zhu  _et al._ [2018] Z. Zhu, G. Sun, W.-L. You, and D.-N. Shi, Fidelity and criticality of a quantum Ising chain with long-range interactions, Phys. Rev. A 98, 023607 (2018).
  * Koziol  _et al._ [2021] J. A. Koziol, A. Langheld, S. C. Kapfer, and K. P. Schmidt, Quantum-critical properties of the long-range transverse-field Ising model from quantum Monte Carlo simulations, Phys. Rev. B 103, 245135 (2021).
  * Haegeman  _et al._ [2011] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Time-Dependent Variational Principle for Quantum Lattices, Phys. Rev. Lett. 107, 070601 (2011).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Unifying time evolution and optimization with matrix product states, Phys. Rev. B 94, 165116 (2016).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011).
  * Browaeys  _et al._ [2016] A. Browaeys, D. Barredo, and T. Lahaye, Experimental investigations of dipole–dipole interactions between a few Rydberg atoms, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 152001 (2016).
  * Tran  _et al._ [2021] M. C. Tran, A. Y. Guo, A. Deshpande, A. Lucas, and A. V. Gorshkov, Optimal State Transfer and Entanglement Generation in Power-Law Interacting Systems, Physical Review X 11, 031016 (2021).
  * Kuwahara and Saito [2020] T. Kuwahara and K. Saito, Strictly Linear Light Cones in Long-Range Interacting Systems of Arbitrary Dimensions, Physical Review X 10, 031010 (2020).
  * Zurek [1985] W. H. Zurek, Cosmological experiments in superfluid helium?, Nature 317, 505 (1985).
  * Zurek [1993] W. H. Zurek, Cosmic strings in laboratory superfluids and the topological remnants of other phase transitions, Acta Phys. Pol. B 24, 1301 (1993).
  * Zurek [1996] W. H. Zurek, Cosmological experiments in condensed matter systems, Phys. Rep. 276, 177 (1996).
  * Zurek and Dorner [2008] W. H. Zurek and U. Dorner, Phase transition in space: how far does a symmetry bend before it breaks?, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 366, 2953 (2008).
  * Damski and Zurek [2009] B. Damski and W. H. Zurek, Quantum phase transition in space in a ferromagnetic spin-1 bose–einstein condensate, New Journal of Physics 11, 063014 (2009).
  * Calabrese and Cardy [2005] P. Calabrese and J. Cardy, Evolution of entanglement entropy in one-dimensional systems, Journal of Statistical Mechanics: Theory and Experiment 2005, P04010 (2005).
  * Calabrese and Cardy [2009] P. Calabrese and J. Cardy, Entanglement entropy and conformal field theory, Journal of Physics A: Mathematical and Theoretical 42, 504005 (2009).
  * Hauke and Tagliacozzo [2013] P. Hauke and L. Tagliacozzo, Spread of correlations in long-range interacting quantum systems, Phys. Rev. Lett. 111, 207202 (2013).
  * Fradkin [2013] E. Fradkin,  _Field theories of condensed matter physics_ (Cambridge University Press, 2013) Chap. 4.2.5.
  * Peskin [1995] M. E. Peskin,  _An introduction to quantum field theory_ (Westview press, 1995).
  * Alonso and De Vincenzo [1997] V. Alonso and S. De Vincenzo, General boundary conditions for a Dirac particle in a box and their non-relativistic limits, Journal of Physics A: Mathematical and General 30, 8573 (1997).
  * Chung  _et al._ [2012] D. J. Chung, L. L. Everett, H. Yoo, and P. Zhou, Gravitational fermion production in inflationary cosmology, Physics Letters B 712, 147 (2012).
  * Gaafar  _et al._ [2019] M. A. Gaafar, T. Baba, M. Eich, and A. Y. Petrov, Front-induced transitions, Nature Photonics 13, 737 (2019).
  * Stepanov [1993] N. S. Stepanov, Waves in nonstationary media, Radiophysics and Quantum Electronics 36, 401 (1993).