#  Supplementary Information

A. De<EMAIL_ADDRESS>Joint Quantum Institute and Joint Center for Quantum
Information and Computer Science, NIST/University of Maryland, College Park,
Maryland 20742  P. Cook<EMAIL_ADDRESS>Department of Physics and Astronomy,
Michigan State University, East Lansing, Michigan 48824  Facility for Rare
Isotope Beams, East Lansing, Michigan 48824  K. Collins  Joint Quantum
Institute and Joint Center for Quantum Information and Computer Science,
NIST/University of Maryland, College Park, Maryland 20742  W. Morong111current
address: AWS Center for Quantum Computing, Pasadena, California 91125, USA.
Work done prior to joining AWS. Joint Quantum Institute and Joint Center for
Quantum Information and Computer Science, NIST/University of Maryland, College
Park, Maryland 20742  D. Paz  Department of Physics and Astronomy, Michigan
State University, East Lansing, Michigan 48824  P. Titum  Joint Quantum
Institute and Joint Center for Quantum Information and Computer Science,
NIST/University of Maryland, College Park, Maryland 20742  Johns Hopkins
University Applied Physics Laboratory, Laurel, Maryland 20723  G. Pagano
Department of Physics and Astronomy, Rice University, Houston, TX, USA  A. V.
Gorshkov  Joint Quantum Institute and Joint Center for Quantum Information and
Computer Science, NIST/University of Maryland, College Park, Maryland 20742
M. Maghrebi  Department of Physics and Astronomy, Michigan State University,
East Lansing, Michigan 48824  C. Monroe  Joint Quantum Institute and Joint
Center for Quantum Information and Computer Science, NIST/University of
Maryland, College Park, Maryland 20742  Department of Electrical and Computer
Engineering, Department of Physics, Duke Quantum Center, Duke University,
Durham, North Carolina 27708, USA  IonQ Inc., College Park, Maryland 20742,
USA

###### Contents

  1. I Dynamical Phase Diagram from Mean-Field Analysis
  2. II LMG model
    1. II.1 Gaussian Fluctuations
    2. II.2 Critical Finite-Size Scaling
    3. II.3 Finite-Size Behavior of the Order Parameter
    4. II.4 Double (Multiple) Quenches
  3. III Power-law Interacting Spin Chain
    1. III.1 Power-Law $J_{ij}$ with Periodic Boundary Conditions
    2. III.2 Power-Law $J_{ij}$ with Open Boundary Conditions
  4. IV Long-range Interacting Spin Chain with Experimental $J_{ij}$
  5. V Experiment vs Theory
    1. V.1 Finite-Size Behavior of the Order Parameter
    2. V.2 Long-Range Character of $J_{ij}$ Matrix
    3. V.3 Correlation Matrix
    4. V.4 Double Quench Switch Times
  6. VI Optimal Collapse Algorithm
    1. VI.1 Uncertainty Estimates via the Hessian

This document provides details of the theoretical analysis (sections I, II,
III and IV), elaborates on and contrasts against the experimental results
(section V), and finally provides the numerical approach used to determine the
critical exponents and their errors (section VI).

## Theoretical Models

In this Supplemental Information, we consider theoretical models with a
generalized long-range interaction:

$H=-\frac{1}{2\mathcal{J}}\sum_{i\neq
j}^{N}J_{ij}\left(\gamma^{x}\sigma_{i}^{x}\sigma_{j}^{x}+\gamma^{y}\sigma_{i}^{y}\sigma_{j}^{y}\right)+B\sum_{i}^{N}\sigma_{i}^{z},$
(1)

where $\sigma_{i}^{\alpha}$ with $\alpha=x,y,z$ are the usual Pauli operators
for spin $i$, and

$\mathcal{J}=\frac{1}{N-1}\sum_{i\neq j}^{N}J_{ij},$ (2)

is a normalization constant known as the Kac factor (or, Kac normalization).
This normalization factor is particularly useful for long-range coupling as it
renders the Hamiltonian extensive.

  1. 1.

In section I, we provide a mean-field analysis to determine the dynamical
phase diagram. We then consider three versions of the model described by eq. 1
ranging from the most idealized theoretical model to one that is most
experimentally relevant.

  2. 2.

In section II, we consider the collective Lipkin-Meshkov-Glick (LMG) model
where $J_{ij}=1$ and the interaction is infinite ranged. We provide a detailed
analysis of the ground-state, the dynamical phase diagram as well as the
critical properties. While an idealized model, it explains the main features
of the experimental results.

  3. 3.

In section III, we add another layer of complexity by considering a power-law
interaction $J_{ij}\sim 1/|i-j|^{p}$ with open and periodic conditions. We
show that, for $0<p<1$, this model exhibits the same dynamical phase diagram
and critical properties.

  4. 4.

Finally, in section IV, we theoretically analyze the long-range spin model
with the experimental values of $J_{ij}$. The latter values fall off
exponentially at large distances (see section V.2), hence departing from the
collective or the truly long-range variant of the previous models.
Nonetheless, we show that the critical properties of the experimental model
are consistent with those of the LMG model and its long-range variants.

Before continuing, we note that the experimental model is antiferromagnetic
while the theoretical models considered here are ferromagnetic; however, we
use a simple trick to relate the two models Islam2013 ; Jurcevic2017 ;
tan2021domain . Given that the Hamiltonian is purely real (using the standard
representation of Pauli matrices) and the fact that the correlation function
$\langle\sigma^{x}_{i}\sigma^{x}_{j}\rangle$ is real, we have

$\langle\sigma^{x}_{i}\sigma^{x}_{j}\rangle_{H,\psi_{0},t}=(\langle\sigma^{x}_{i}\sigma^{x}_{j}\rangle_{H,\psi_{0},t})^{*}=\langle\sigma^{x}_{i}\sigma^{x}_{j}\rangle_{-H,\psi_{0},t},$
(3)

where we have made explicit the dependence on the Hamiltonian $H$, time $t$
and the initial state $\psi_{0}$. We have furthermore assumed that the initial
state is also real; for example, a fully polarized state along the negative
$z$ direction is real. By flipping the overall sign of the Hamiltonian, the
antiferromagnetic coupling is then mapped to a ferromagnetic one.

##  I Dynamical Phase Diagram from Mean-Field Analysis

In this section, we obtain the dynamical phase diagram using a mean-field
analysis for the generalized model in eq. 1. Starting from the initial state
$\left|\downarrow\downarrow\cdots\right\rangle$, we suddenly turn on an Ising-
type Hamiltonian by setting $\gamma^{x}=1$ and $\gamma^{y}=0$. In a mean-field
analysis, the Hamiltonian becomes

$H=-m\sigma^{x}+B\sigma^{z}=\sqrt{m^{2}+B^{2}}\,\,{\mathbf{n}}\cdot{\boldsymbol{\sigma}},$
(4)

where $m\equiv\langle\sigma^{x}\rangle$. We have dropped the site index, and
defined the unit vector $\mathbf{n}=(n_{x},0,n_{z})$ parallel to $(-m,0,B)$.
At long times after the sudden quench, the initial state dephases in the new
basis defined by the operator ${\mathbf{n}}\cdot{\boldsymbol{\sigma}}$, and
the time-averaged density matrix becomes

$\lim_{t\to\infty}|\psi\rangle_{t}\langle\psi|=\sin^{2}\frac{\theta}{2}|{\mathbf{n}},+\rangle\langle{\mathbf{n}},+|+\cos^{2}\frac{\theta}{2}|{\mathbf{n}},-\rangle\langle{\mathbf{n}},-|,$
(5)

where $\mathbf{n}\cdot\mathbf{z}=\cos\theta$. Using the identity
$m=\langle\sigma^{x}\rangle$, we then find the self-consistent equation

$m=-\cos\theta\sin\theta=\frac{m}{\sqrt{m^{2}+B^{2}}}\frac{B}{\sqrt{m^{2}+B^{2}}},$
(6)

which in turn yields the time-averaged order parameter as

$m=\sqrt{B\left(1-B\right)}.$ (7)

Furthermore, this equation predicts a dynamical phase transition at $B_{c}=1$.
This mean-field analysis is exact for the infinite-range LMG model.
Furthermore, with the aid of spin-wave analysis, we show that it is also exact
for the long-range model with $0<p<1$ (see section III), and additionally
provides an excellent estimate for the phase transition corresponding to the
experimental $J_{ij}$ (see section IV).

##  II LMG model

The LMG model is characterized by the infinite-range interaction with
$J_{ij}\equiv 1$. This model can be reformulated in terms of total spin
operators, with the Hamiltonian (up to a constant)

$H_{LMG}=-\frac{2}{N}\left(\gamma^{x}S_{x}^{2}+\gamma^{y}S_{y}^{2}\right)+2BS_{z},$
(8)

where $S_{\alpha}=\frac{1}{2}\sum_{i}\sigma^{\alpha}_{i}$. At zero
temperature, this model exhibits a quantum phase transition at $B_{gs,c}=1$
suzuki2012quantum . Interestingly, the ground state critical point coincides
exactly with the dynamical critical point obtained in the previous section.
However, there are important distinctions: the order parameter’s dependence on
$B$ is different. In fact, the ground state is perfectly ordered at $B=0$
while the corresponding order parameter long after the sudden quench vanishes
when $B=0$. More importantly, the critical behaviour at the respective phase
transitions is different: while the ground state exhibits a quantum phase
transition at $B=1$, the dynamical critical point shows distinct (and,
effectively thermal) critical behaviour (see section II.2). Furthermore, we
show that for a double quench, the critical behaviour is neither quantum nor
thermal, and is genuinely non-equilibrium (see section II.4).

### Beyond Mean Field

Next, we consider fluctuations at the phase transition. While the mean field
solution in eq. 7 exactly describes the ordered phase (where $B<B_{c}=1$) in
the thermodynamic limit, fluctuations and finite-size scaling become important
at or near the critical point for a finite system. To characterize
fluctuations in the disordered phase as well as critical fluctuations near the
critical point, one can use the Holstein-Primakoff transformation to map the
spin operators of the LMG model to bosonic creation and annihilation operators
as ($S_{\pm}=S_{x}+iS_{y}$)

$S_{-}=S_{+}^{\dagger}=\sqrt{N-a^{\dagger}a}a\qquad S_{z}=-N/2+a^{\dagger}a.$
(9)

To the lowest order in $1/N$, we can identify

$S_{x}\approx\sqrt{\frac{N}{2}}x,\qquad S_{y}\approx\sqrt{\frac{N}{2}}p.$ (10)

Assuming that $\gamma^{x}>\gamma^{y}$, the dominant Ising interaction is along
the $x$ direction. In this case, the LMG Hamiltonian can be simply written as
that of a harmonic oscillator Titum20 ,

$H=\frac{1}{2m}p^{2}+\frac{1}{2}m\Omega^{2}x^{2}+\frac{1}{4N}ux^{4}+\cdots,$
(11)

where we have defined $m^{-1}=2\left(B-\gamma^{y}\right)$,
$\Omega^{2}=4\left(B-\gamma^{x}\right)\left(B-\gamma^{y}\right)$, and
$u=2\gamma^{x}$. The dots in the above equation represent higher-order terms
in $1/N$ or higher powers of $x$ and $p$ which are irrelevant at or away from
the critical point in the disordered phase.

###  II.1 Gaussian Fluctuations

We begin by considering eq. 11 at the quadratic level in the disordered phase,
before taking into account the interaction term ($\propto 1/N$). We are
interested in the dynamics upon the sudden quench
$(\gamma^{x,y}_{0},B_{0})\to(\gamma^{x,y},B)$. Equivalently, we can consider a
quench $(m_{0},\Omega_{0})\longrightarrow(m,\Omega)$ of the harmonic
oscillator, where it is easy to solve for the time-averaged fluctuations:

$\frac{1}{N}\overline{\langle
S_{x}^{2}\rangle}=\frac{1}{2}\overline{\left\langle
x^{2}\right\rangle}=\frac{m_{0}\Omega_{0}}{8m^{2}\Omega^{2}}\left(1+\frac{m^{2}\Omega^{2}}{m_{0}\Omega_{0}^{2}}\right)\,,$
(12)

where the overline indicates time averaging. Notice that fluctuations diverge
where $\Omega=0$, which in turn sets $(B-\gamma^{x})(B-\gamma^{y})=0$. This
condition coincides with the mean-field prediction in eq. 7 when
$\gamma^{x}=1$ and $\gamma^{y}=0$; this is expected as the LMG model is
infinite ranged. Furthermore, we can characterize the critical fluctuations
beyond the mean-field prediction. We first observe that fluctuations scale as
$\overline{\langle x^{2}\rangle}\sim 1/\Omega^{2}$ near the critical point.
Comparing against the expression for a harmonic oscillator at high
temperatures

$\frac{1}{2}m\Omega^{2}\left\langle x^{2}\right\rangle\stackrel{{\scriptstyle
T\gg\Omega}}{{\approx}}\frac{1}{2}T,$

we can identify an effective temperature Titum20

$T_{\mathrm{eff}}=\frac{m_{0}\Omega_{0}}{4m}=\frac{B}{2}\sqrt{\frac{B_{0}-\gamma^{x}_{0}}{B_{0}-\gamma^{y}_{0}}}.$
(13)

Specifically, in a quench starting from a product state corresponding to the
ground state in the absence of the Ising interaction ($\gamma^{x,y}_{0}=0$) to
the critical point of the infinite-range Ising Hamiltonian ($\gamma^{x}=1$,
$\gamma^{y}=0$, and $B=1$), we find a constant effective temperature, $T_{\rm
eff}=1/2$. We stress that such notion of the effective temperature does not
mean that the system is in equilibrium; rather, it means that low-frequency
modes (or, more precisely, the soft modes) have become effectively thermal.
Indeed, one finds thermal critical behaviour and scaling in this type of
quench Titum20 ; see also fig. S1.

###  II.2 Critical Finite-Size Scaling

Fluctuations diverge as we approach the critical point $B\to 1$ either in the
ground state or in a quench to the critical point. In a finite-size system,
the fluctuations do not strictly diverge but rather scale with $N$ in a
nontrivial fashion:

$\frac{1}{N}\langle S_{x}^{2}\rangle\sim N^{\alpha}.$ (14)

Notice that in the disordered phase, we have $\alpha=0$ indicating noncritical
fluctuations while in the ordered phases $\alpha=1$ due to ordering $\langle
S_{x}\rangle\sim N$. Only at the critical point, $\alpha$ becomes a nontrivial
critical exponent. For the quantum phase transition in the ground state, this
exponent becomes $\alpha=1/3$; see fig. S1(a). This should be contrasted with
a thermal phase transition where this exponent becomes $\alpha=1/2$. In a
quench to the critical point, the system evolves and approaches a stationary
state that is not an equilibrium state; however, the critical exponent takes
the same value $\alpha=1/2$ as thermal equilibrium; see fig. S1(b). One can
view this as a consequence of the fact that the effective temperature takes a
finite nonzero value (see eq. 13). Here, we derive this critical scaling using
a simple analysis. Our approach is based on a non-equilibrium field theory
which systematically leads to a semiclassical analysis similar to the
truncated Wigner approximation (see Ref. Titum20  for more details). To this
end, we start from the partition function

$Z=\int Dx_{\pm}W_{0}e^{i(S[x_{+}]-S[x_{-}])}\,,$ (15)

where $W_{0}$ is the Wigner function corresponding to the initial state, and
$[x_{\pm}(t)]$ denote forward and backward trajectories, respectively. These
trajectories are weighted by the phase factors $\exp(\pm iS)$ where
$S=\frac{1}{2}\int dt[m\dot{x}^{2}-m\Omega^{2}x^{2}-(u/2N)x^{4}]$ is the
corresponding action. Introducing the (Keldysh) variables $x_{c/q}=(x_{+}\pm
x_{-})/\sqrt{2}$, the total (Keldysh) action becomes

$S_{K}=S[x_{+}]-S[x_{-}]=-\int
dt[mx_{q}\ddot{x}_{c}+rx_{q}x_{c}+(u/2N)(x_{c}^{3}x_{q}+x_{c}x_{q}^{3})]\,,$
(16)

with $r=m\Omega^{2}$. We take the initial state as the ground state of the
Hamiltonian in the absence of the Ising interaction. The Wigner function then
becomes a Gaussian function of phase space variables as

$W_{0}(x_{0},p_{0})\sim e^{-x_{0}^{2}-p_{0}^{2}}\,.$ (17)

Casting the Wigner function in terms of the configuration variables, the total
partition function takes the form

$Z=\int Dx_{c/q}e^{-x_{c0}^{2}/2l_{0}^{2}-\dot{x}_{c0}^{2}/2v_{0}^{2}}e^{-\int
dt[mx_{q}\ddot{x}_{c}+rx_{q}x_{c}+(u_{c}/2N)x_{c}^{3}x_{q}+(u_{q}/2N)(x_{c}x_{q}^{3})]}\,,$
(18)

where we have introduced the parameters $l_{0}$ and $v_{0}$ for later
convenience; for the coherent initial state assumed above, we have $l_{0}=1$
and $v_{0}=1/m_{0}$. We have also introduced the coefficients $u_{c/q}$ to
distinguish the two interaction vertices; at the microscopic level, we have
$u_{c}=u_{q}=u$. As pointed out in Ref. Titum20 , to probe the behaviour at
long times, we can rescale time $t\to t/\lambda$ and identify how various
terms scale. We find that that under this rescaling $l_{0}\to 0$ and $u_{q}\to
0$. The first condition imposes the boundary condition $x_{c0}=0$ given the
Gaussian distribution $\exp(-x_{c0}^{2}/2l_{0}^{2})$. The latter condition
simply means that the ‘quantum’ vertex proportional to $u_{q}$ is less
relevant than the classical vertex proportional to $u_{c}$, a condition that
typically emerges in the semiclassical limit. Dropping the quantum vertex, we
can integrate over the field $x_{q}$ to find a delta function that imposes the
semiclassical equation of motion together with a Gaussian distribution of the
initial state (also changing the coordinate back to $x$ via $x_{c}=\sqrt{2}x$
and restoring $u_{c}\to u$):

$Z=\int
Dx\delta(x_{0})e^{-\dot{x}_{0}^{2}/v_{0}^{2}}\delta\left(m\ddot{x}+rx+\frac{u}{N}x^{3}\right)\,.$
(19)

This equation admits a simple interpretation: A ‘particle’ whose position is
given by $x$ moves under a nonlinear equation of motion with the initial
conditions where $x_{0}=0$ and $\dot{x}_{0}$ is drawn from a Gaussian
distribution. This equation is simply the truncated Wigner approximation,
which evolves the wavefunction with the classical equation of motion, but
captures the (quantum) fluctuations of the initial state; the only difference
is that, using a scaling argument, we have set $x_{0}=0$.

To characterize the critical behaviour, we can study the scaling behaviour of
the above equation. We remark that, at the critical point $r=0$, the partition
function becomes scale invariant under the rescaling

$t\to t/\lambda,\quad x\to x/\lambda,\quad N\to N/\lambda^{4}\,.$ (20)

This immediately implies that fluctuations $\langle x^{2}\rangle$ as a
function of $t$ and $N$ both take the scaling form

$\langle x(t)^{2}\rangle=N^{\alpha}f(t/N^{\zeta})\,,$ (21)

with the critical exponents $\alpha=1/2$ and the $\zeta=1/4$; the latter
exponent identifies a dynamical critical exponent. This is indeed consistent
with the exact numerical calculation in fig. S1(b).

(a)  (b)  (c)

Figure S1: Fluctuations scaling with system size at the critical point (a) in
the ground state, (b) after a single quench, and (c) after a double quench.
(b) A single quench from the disordered phase to the critical point
${\left(B,\gamma^{x},\gamma^{y}\right)=(1,1,0)}$ gives rise the critical
exponents ${\left(\alpha_{1},\zeta_{1}\right)=\left(1/2,1/4\right)}$. (c) A
second quench to the other critical point
${\left(B,\gamma^{x},\gamma^{y}\right)=(1,0,1)}$ is performed when the
fluctuations in panel (b) reach their peak. The coordinate $t_{2}$ denotes the
time since the second quench. The double quench gives rise to the critical
exponents ${\left(\alpha_{2},\zeta_{2}\right)=\left(3/4,1/8\right)}$.

###  II.3 Finite-Size Behavior of the Order Parameter

In eq. 7, we have presented the mean-field prediction for the time-averaged
order parameter in the ordered phase, $B<1$, and in the thermodynamic limit
$N\to\infty$. In the previous section, we have further investigated the
critical fluctuations for finite system sizes exactly at the critical point
$B=B_{c}=1$. In this section, we study finite-size corrections in the vicinity
of the critical point $B_{c}$. This is crucial to identify the dynamical phase
transition in a finite-size system. To this end, we consider the equation of
motion dictated by eq. 19:

$\ddot{x}+rx+ux^{3}=0.$ (22)

For notational convenience, we have set $m=1$ and absorbed the factor
proportional to $1/N$ into the definition of the interaction parameter $u$.
For the purpose of this section, it suffices that

$r\propto B-1,\qquad u\propto 1/N\,.$ (23)

In the quench starting from a disorded state, the initial conditions are set
by $x(t=0)=0$ and $\dot{x}(0)=v$ where $v$ is drawn from a Gaussian
distribution $\exp(-Dv^{2})$; the parameter $D$ is set by the initial state,
but we shall treat it as a phenomenological parameter. Solving the above
equation with the initial conditions $x(t=0)=0$ and $\dot{x}(t=0)=v$, we find

$x_{v}(t)=\mathrm{sgn}(v)\,{\mathrm{sn}}\left(t\sqrt{\frac{r+\sqrt{r^{2}+2uv^{2}}}{2}},\frac{r-\sqrt{r^{2}+2uv^{2}}}{r+\sqrt{r^{2}+2uv^{2}}}\right)\,,$
(24)

with ${\rm sn}(u,m)$ the Jacobi elliptic function. The above equation exhibits
persistent oscillations. We are instead interested in the time average of
$\overline{x(t)^{2}}$ which is given by

$\overline{x_{v}^{2}}\equiv\lim_{t\to\infty}\overline{x_{v}(t)^{2}}=\frac{\sqrt{r^{2}+2uv^{2}}-r}{2u}\,.$
(25)

Finally, the Gaussian integral over $v$ (i.e., the integral
$\langle\cdot\rangle_{v}\equiv\int dv\exp(-Dv^{2})\cdot/\int dv\exp(-Dv^{2})$)
yields

$\left\langle\overline{x_{v}^{2}}\right\rangle_{v}=-\frac{r}{2u}+\frac{U(-\frac{1}{2},0,Dr^{2}/(2u))}{\sqrt{2Du}}\,,$
(26)

where $U$ is Tricomi’s confluent hypergeometric function. Making the
dependence of $u$ on $N$ explicit by $u\to u/N$, we find

$\overline{{m}^{2}}\equiv\frac{1}{N^{2}}\overline{\langle
S_{x}(t)^{2}\rangle}=\frac{1}{N}\left\langle\overline{x_{v}^{2}}\right\rangle_{v}=-\frac{r}{2u}+\frac{1}{\sqrt{N}}\frac{U\left(-\frac{1}{2},0,DNr^{2}/(2u)\right)}{\sqrt{2Du}}\,.$
(27)

Deep in the ordered phase, we must have $\overline{m^{2}}=m^{2}$ where the
order parameter $m$ is computed from mean field (see eq. 7). This is indeed
the case: the above equation captures the initial linear dependence of the
order parameter on $r$ as one enters the ordered phase, that is,
$\overline{m^{2}}\sim|r|\propto(1-B)$ when $r<0$, while the exact (mean-field)
solution within the ordered phase is given by $m^{2}=B(1-B)$. An expression
that properly interpolates between the two extremes can be obtained by
multiplying the right hand side of the above equation by a factor of $B$. We
also fix the coefficient $u=1/2$ to match eq. 7 deep in the ordered phase
where $r<0$. Making explicit the dependence on $B$ and $B_{c}$, we find

$\overline{m^{2}}=\frac{B}{B_{c}}\left[\left(1-\frac{B}{B_{c}}\right)+\frac{1}{\sqrt{N}}\frac{U\left(-\frac{1}{2},0,DN(1-\frac{B}{B_{c}})^{2}\right)}{\sqrt{D}}\right]\,.$
(28)

In a numerical calculation of the critical point, the above fit only contains
two free parameters: the constant $D$ and the critical point $B_{c}$.

In the experiment, it is more convenient to study the maximum of $\langle
S_{x}^{2}\rangle/N^{2}$ in time rather than its stationary value. The above
analysis remains valid but with an an overall amplitude, which we shall treat
as a free parameter in eq. 28. In fig. S2, we show that that above function
(up to an overall amplitude $\cal A$) is an excellent fit to the order
parameter obtained from the exact numerical simulation, and the extracted
value of the critical point is almost identical to the exact value of
$B_{c}=1$.

Figure S2: The order parameter for the LMG model with $N=50$ extracted from
the peak of $\langle S_{x}^{2}(t)\rangle$. We also depict the fit function
from eq. 28 with an overall factor $\mathcal{A}=16.15(9)$, as well as the fit
parameters $D=0.34(1)$ and $B_{c}=0.976(4)$ in excellent agreement with the
exact critical point, $B_{c}=1$.

###  II.4 Double (Multiple) Quenches

In this section, we analyze the scaling behaviour of the system as the result
of a double quench starting from a disordered initial state to the critical
point of the LMG model with $\gamma^{x}_{0}=1$ and $\gamma^{y}_{0}=0$, and,
with a delay, a subsequent quench to the critical point of the LMG model with
$\gamma^{x}_{0}=0$ and $\gamma^{y}_{0}=1$; See Fig. 1(b) of the main text.
Both critical points are defined by $B_{c}=1$. Notice that the final
Hamiltonian has the Ising interaction along the $y$ direction only, and the
dominant fluctuations correspond to the those along the $y$ direction. To
parallel the notation in the previous sections, we use the Holstein-Primakoff
expansion but now identify, to the lowest order in $1/N$,

$\sqrt{\frac{2}{N}}S_{x}\approx p_{2}\,,\qquad\sqrt{\frac{2}{N}}S_{y}\approx
x_{2}\,,$ (29)

as the dynamical variables upon a second quench, while we reserve $x_{1}$ and
$p_{1}$ defined similarly to $x$ and $p$ in eq. 10, respectively, as the
dynamical variables following the first but before the second quench. At the
time of the second quench ($t_{2}=0$), we have $x_{2}=p_{1}$ and $p_{2}=x_{1}$
due to continuity. The Hamiltonian at each stage is given by
$H_{i}=\frac{1}{2m_{i}}p_{i}^{2}+\frac{1}{2}m_{i}\Omega_{i}^{2}x_{i}^{2}$
where $i=0,1,2$ correspond to pre-quench, first-quench, and second-quench
variables; we also identify $x_{0}=x_{1}$ and $p_{0}=p_{1}$. To characterize
fluctuations at the quadratic order, we shall consider a finite distance from
the critical point at all stages of the dynamics (but approaching criticality
when considering finite-size scaling). The frequencies $\Omega_{i}$ are then
given by

$\Omega_{i}^{2}=4\left(B_{i}-J\gamma^{x}_{i}\right)\left(B_{i}-J\gamma^{y}_{i}\right)$
(30)

and similarly for $m_{i}$; the exact values of $m_{i}$ are unimportant and we
just remark that they remain finite even as we approach the critical point.

Next, we characterize the dominant fluctuations ($\langle x_{2}^{2}\rangle$)
at late times from those of $\langle x_{1}^{2}\rangle$ and $\langle
p_{1}^{2}\rangle$ long after the first quench but before the second quench.
Finally taking the long-time average, we find

$\overline{\left\langle x_{2}^{2}\right\rangle}=\frac{\overline{\left\langle
p_{1}^{2}\right\rangle}}{2}+\frac{\overline{\left\langle
x_{1}^{2}\right\rangle}}{2m_{2}^{2}\Omega_{2}^{2}}=\frac{m_{0}^{2}\Omega_{0}^{2}}{16}+\frac{m_{1}^{2}\Omega_{1}^{2}}{16m_{0}^{2}\Omega_{0}^{2}}+\frac{1}{8m_{0}\Omega_{0}m_{2}^{2}\Omega_{2}^{2}}+\frac{m_{0}\Omega_{0}}{8m_{1}^{2}\Omega_{1}^{2}m_{2}^{2}\Omega_{2}^{2}}.$
(31)

where the frequencies $\Omega_{i}$ are defined above. For two successive
quenches to the vicinity of critical points,
$\Omega_{2}\ll\Omega_{1}\ll\Omega_{0}$, the above expression reduces to

$\frac{1}{2}\overline{\left\langle
x_{2}^{2}\right\rangle}\sim\frac{m_{0}\Omega_{0}}{16m_{1}^{2}\Omega_{1}^{2}m_{2}^{2}\Omega_{2}^{2}}.$
(32)

Using the equipartition theorem $\frac{1}{2}m_{2}\Omega_{2}^{2}\langle
x_{2}^{2}\rangle\sim\frac{1}{2}T$, we can now identify an effective
temperature upon the second quench (denoted by subscript 2) as

$T_{\mathrm{eff},2}=\frac{m_{0}\Omega_{0}}{8m_{2}m_{1}^{2}\Omega_{1}^{2}}.$
(33)

Notice, however, that the effective temperature diverges as we tune the first
quench to the critical point, $\Omega_{1}\to 0$. Such divergence hints at a
critical behaviour that is strongly athermal (see also Ref. Titum20 ).

To characterize the critical behaviour, we can set up a scaling analysis using
semiclassical dynamics. Using the continuity of these functions at $t_{2}=0$
where $x_{2}=p_{1}$ and $\dot{x}_{2}\sim p_{2}=x_{1}$, the initial conditions
for $x_{2}$ and its time derivative are given by a distribution function; we
can alternatively characterize the same information via the expectation values

$\langle x_{2}^{2}\rangle\sim N^{0},\quad\langle\dot{x}_{2}^{2}\rangle\sim
N^{\alpha_{1}},\qquad\mbox{at}\,\,t_{2}=0\,.$ (34)

For a uniform notation, we have identified $\alpha_{1}=\alpha$ as the critical
exponent emerging in the first quench. To this end, a rather similar scaling
analysis gives the partition function

$Z=\int
Dx_{2}\delta[x_{2}(0)]P\left[(\dot{x}_{2}(0))^{2}/N^{\alpha_{1}}\right]\delta\left(m_{2}\ddot{x}_{2}+r_{2}x_{2}+\frac{u_{2}}{N}x_{2}^{3}\right)\,,$
(35)

where $r_{2}=m_{2}\Omega_{2}^{2}$ and $u_{2}=2\gamma^{y}_{2}$; we have also
defined $x_{2}(0)\equiv x_{2}(t_{2}=0)$ for ease of notation (similarly for
$\dot{x}_{2}(0)$). To arrive at this equation, we have used the scaling limit
where the quantum vertex can be ignored and the fluctuations of $x_{2}$ right
after the second quench ($t_{2}=0$) becomes a delta function. On the other
hand, we keep track of the initial fluctuations of $\dot{x}_{2}$ and its
scaling with system size $N$ via the the distribution function $P$; for our
purposes, only the scaling behaviour, and not the exact form of the function
$P$, is important. We stress that the scaling variable in the argument of the
function $P$ is dictated by scaling behaviour in eq. 34 at the onset of the
second quench.

Next, we observe that the partition function is scale-invariant (at $r_{2}=0$)
upon the transformation

$t\to t/\lambda,\quad x\to x/\lambda^{\frac{1+\alpha_{1}}{1-\alpha_{1}}},\quad
N\to N/\lambda^{\frac{4}{1-\alpha_{1}}}\,.$ (36)

We can then immediately write a scaling relation for fluctuations

$\langle x_{2}(t_{2})^{2}\rangle=N^{\alpha_{2}}g(t_{2}/N^{\zeta_{2}})\,,$ (37)

with $g$ a scaling function and the critical exponents identified as

$\alpha_{2}=\frac{1+\alpha_{1}}{2},\qquad\zeta_{2}=\frac{1-\alpha_{1}}{2}\,.$
(38)

Using $\alpha_{1}=1/2$ in the first quench, we find $\alpha_{2}=\frac{3}{4}$
and $\zeta_{2}=\frac{1}{8}$. These predictions are in excellent agreement the
exact numerical simulation shown in fig. S1(c).

At this point, we can also determine the scaling of the effective temperature.
Note that eq. 35 remains scale-invariant even for nonzero $r_{2}$ upon scaling
$r_{2}\to\lambda^{2}r_{2}$, which then determines the scaling behaviour of
this variable as $r_{2}\sim N^{-2\zeta_{2}}$. Now recalling the definition of
the effective temperature $\frac{1}{2}\Omega_{2}^{2}\overline{\langle
x_{2}^{2}\rangle}\sim\frac{1}{2}T_{\rm eff,2}$ together with the facts
$r_{2}\propto\Omega_{2}^{2}\sim N^{-2\zeta_{2}}$ and $\langle
x_{2}^{2}\rangle\sim N^{\alpha_{2}}$, we find

$T_{\rm eff,2}=N^{\alpha_{2}-2\zeta_{2}}\,.$ (39)

The same relation holds for the first quench by changing all the subscripts to
1. While in the latter case, $\alpha_{1}-2\zeta_{1}=0$ indicating a constant
effective temperature, for a double quench $\alpha_{2}-2\zeta_{2}=1/2$
highlighting a nontrivial scaling of the effective temperature and a genuinely
nonthermal critical behaviour.

We also remark that more generally for $k+1$ consecutive critical quenches
$\gamma^{x}\to\gamma^{y}\to\gamma^{x}\to\cdots$ (showing only the nonzero
anisotropy parameter) starting from a disordered state, we find the non-
equilibrium critical exponents:

$\alpha_{k+1}=\frac{1+\alpha_{k}}{2},\qquad\zeta_{k+1}=\frac{1-\alpha_{k}}{4}\,,$
(40)

and the effective temperature scaling

$T_{\mathrm{eff},k}\sim N^{\alpha_{k}-2\zeta_{k}}\,,$ (41)

with $\alpha_{0}=0$. One can solve the above iterative equations to find the
exponents

$\alpha_{k}=1-2^{-k},\qquad\zeta_{k}=2^{-k-1}\,.$ (42)

The critical exponents for single and double quench are recovered by setting
$k=1$ and 2, respectively. We can also determine the critical scaling of the
temperature from eq. 41 as

$T_{\mathrm{eff},k}\sim N^{1-\frac{1}{2^{k-1}}}\,.$ (43)

##  III Power-law Interacting Spin Chain

In this section, we extend our analysis to include the long-range model in eq.
1 where the interaction $J_{ij}$ falls off as a power-law with the distance,
$J_{ij}\sim 1/|i-j|^{p}$. We show, using a spin-wave analysis, that the
universal behaviour discussed in the text is governed by a collective mode as
long as $0<p<1$, hence the same universality class as the LMG model discussed
in the previous section. Furthermore, we show that the critical point is
consistent with the mean-field (or, the LMG model’s) prediction, i.e.,
$B_{c}=1$, in the thermodynamic limit.

We begin by applying the Holstein-Primakoff transformation to eq. 1 at the
level of individual spins

$\displaystyle\sigma^{-}_{i}=(\sigma_{i}^{+})^{\dagger}=\sqrt{1-a^{\dagger}_{i}a_{i}}a_{i}\approx
a_{i},\qquad\sigma^{z}_{i}=1-2a^{\dagger}_{i}a_{i}\,,$ (44)

where we have neglected the non-linear terms. This approximation is valid as
long as we remain in the disordered phase where the total spin is fully
polarized along the $z$-direction, $\langle S_{z}\rangle=N/2$. Excitations of
the bosonic degrees of freedom amount to fluctuations beyond this state. Under
this approximation, the bosonic form of the Hamiltonian is given by

$H=-\frac{1}{2\mathcal{J}}\sum_{ij}J_{ij}(a_{i}+a^{\dagger}_{i})(a_{j}+a^{\dagger}_{j})+2B\sum_{i}a^{\dagger}_{i}a_{i}\,,$
(45)

where we have neglected an unimportant constant, and have only included the
quadratic terms in the Hamiltonian. This Hamiltonian, being quadratic in the
bosonic operators, can be diagonalized via a Bogoliubov transformation. The
form of this transformation depends on the form of $J_{ij}$. We first consider
power-law interactions with periodic boundary conditions where simple
analytical expressions can be obtained, and then extend our results to open
boundary conditions. In the next section, we apply the same analysis to the
experimental interaction matrix and contrast the results against those with a
purely power-law interaction.

###  III.1 Power-Law $J_{ij}$ with Periodic Boundary Conditions

For illustrative purposes, we first consider the case where $J_{ij}$ is
described by a power-law subject to periodic boundary conditions:

$J_{ij}=\frac{J}{r_{ij}^{p}},\qquad\mbox{with}\quad
r_{ij}=\text{min}(|i-j|,N-|i-j|)\,,$ (46)

and $J_{ij}=0$ when $i=j$. We can bring the Hamiltonian into a convenient form
by going to the Fourier basis,

$a_{j}=\frac{1}{\sqrt{N}}\sum_{k}e^{ijk}a_{k}\,,$ (47)

where $k=2\pi n/N$, $n\in\left\\{0,\ldots,N-1\right\\}$. The Hamiltonian now
becomes diagonal in $k$,

$H=\sum_{k}\begin{pmatrix}a^{\dagger}_{k}\\\
a_{-k}\end{pmatrix}^{T}\begin{pmatrix}-\frac{\tilde{J}_{k}}{2\mathcal{J}}+B&-\frac{\tilde{J}_{k}}{2\mathcal{J}}\\\
-\frac{\tilde{J}_{k}}{2\mathcal{J}}&-\frac{\tilde{J}_{k}}{2\mathcal{J}}+B\end{pmatrix}\begin{pmatrix}a_{k}\\\
a^{\dagger}_{-k}\end{pmatrix}\,,$ (48)

but remains to be diagonalized in terms of the bosonic operators. Here, we
have defined

$\tilde{J}_{k}=\sum_{r=-(N-1)/2}^{(N-1)/2}J_{r}e^{-irk}=2J\sum_{r=1}^{(N-1)/2}\frac{1}{r^{p}}\cos(rk)\,,$
(49)

and assumed that $N$ is odd for simplicity. Next, we employ the Bogoliubov
transformation,

$a_{k}=\cosh(\theta_{k})b_{k}+\sinh(\theta_{k})b^{\dagger}_{-k}\,,$ (50)

where $\theta_{k}=\theta_{-k}$, and the coefficients have been chosen to
ensure that the new operators $b_{k}$ satisfy the bosonic commutation relation
$[b_{k},b^{\dagger}_{k^{\prime}}]=\delta_{kk^{\prime}}$. To determine
$\theta_{k}$, we demand that this transformation diagonalizes eq. 48 as
$H=\sum_{k}\omega_{k}b^{\dagger}_{k}b_{k}$, which leads to the relation

$\tanh(2\theta_{k})=\frac{\tilde{J}_{k}/2\mathcal{J}}{B-\tilde{J}_{k}/2\mathcal{J}}\,,$
(51)

determining $\theta_{k}$, as well as the mode frequencies

$\omega_{k}=\sqrt{B(B-\tilde{J}_{k}/\mathcal{J})}\,.$ (52)

With $\tilde{J}_{0}/\mathcal{J}=(N-1)/N$, we have $\omega_{0}=\sqrt{B(B-1)}$
in the thermodynamic $N\to\infty$ limit. Therefore, the collective mode
corresponding to $k=0$ becomes gapless (or, softens) at the phase transition,
and should be identified as the soft mode of the phase transition. For $B<1$,
the dispersion becomes complex valued and the Holstein-Primakoff
transformation is no longer valid. Most importantly, we find that all the
other modes beside the $k=0$ mode remain gapped for long-range interactions
with $p<1$. This is because $\tilde{J}_{k}/{\cal J}<1$ for any $k>0$,
therefore $\omega_{k}>0$ at or away from the phase transition within the
normal phase, $B\geq 1$. Indeed, one can show that the gap between the $k=0$
mode and the next mode corresponding to $k=2\pi/N$ remains finite for $p<1$
even in the thermodynamic limit $N\to\infty$ Defenu2021 :

$\lim_{N\to\infty}\omega_{k=2\pi/N}-\omega_{k=0}>0,$ (53)

at or beyond the critical point. This implies that the critical behaviour of
long-range interactions with $p<1$ is in the same universality class of the
LMG model.

###  III.2 Power-Law $J_{ij}$ with Open Boundary Conditions

It is generally expected that the bulk properties of a large system would be
insensitive to boundary conditions at least for short-range interactions. The
role of boundaries are, however, more pronounced in the presence of long-range
interactions. Specifically, for truly long-range models with $p<1$, the
boundary conditions can even change the qualitative behaviour Hauke2013 . In
this section, we show that open boundary conditions do not alter our
conclusions for a long-range interacting chain with periodic boundary
conditions. To perform a spin-wave analysis for $J_{ij}$ (again taking
$J_{ii}=0$) without translation invariance, we cannot use the Fourier basis
and instead carry out the Bogoliubov transformation directly in real space. We
rewrite the Hamiltonian in the more convenient form

$H=\sum_{i,j}\begin{pmatrix}a_{i}\\\
a^{\dagger}_{i}\end{pmatrix}^{\dagger}\begin{pmatrix}A_{ij}&C_{ij}\\\
C_{ij}&A_{ij}\end{pmatrix}\begin{pmatrix}a_{j}\\\
a^{\dagger}_{j}\end{pmatrix}\,,$ (54)

where $A_{ij}=-J_{ij}/2\mathcal{J}+B\delta_{ij}$,
$C_{ij}=-J_{ij}/2\mathcal{J}$. Using the fact that the matrices $\mathbf{A}$
and $\mathbf{C}$ commute (the bold notation denotes matrix or vector objects),
we can directly determine the Bogoliubov transformation which brings the above
Hamiltonian into a diagonal form $H=\sum_{k}\Lambda_{k}b^{\dagger}_{k}b_{k}$
colpa_diagonalization_1978 . Specifically, the mode energies $\Lambda_{k}$ are
given by the square root of the eigenvalues of the matrix

$\mathbf{M}=(\mathbf{A}+\mathbf{C})(\mathbf{A}-\mathbf{C})=\mathbf{A}^{2}-\mathbf{C}^{2}\,.$
(55)

The eigenvalues can be determined numerically for an arbitrary interaction
matrix $J_{ij}$.

In fig. S3, we consider $J_{ij}$ matrix describing a power-law interacting
chain with open boundary conditions and plot the mode frequencies
$\Lambda_{0}^{2}$ and $\Lambda_{1}^{2}$ for $N=50$ and $p=0.9$; the value of
$p$ is chosen the mirror the experimental model although the latter features
an exponential decay at large distances (see the next section). $\Lambda_{0}$
corresponds to the soft mode which, by definition, becomes gapless at the
phase transition, $B\simeq 1$. (We recover $B_{c}=1$ exactly in the
thermodynamic limit.) In addition, we see that $\Lambda_{1}$ remains gapped
even at the critical point, and therefore does not contribute the critical
behaviour. Again, $\Lambda_{0}^{2}$ becomes negative for $B<1$, and our
bosonic mapping is no longer valid.

(a)  Figure S3: The two lowest energy eigenvalues $\Lambda_{0}$ and
$\Lambda_{1}$ corresponding to a power-law interaction with open boundary
conditions, $p=.9$ for a system size of $N=50$. While the lowest-frequency
mode softens (i.e., $\Lambda_{0}\to 0$) at the critical point, all the other
modes remain gapped.

##  IV Long-range Interacting Spin Chain with Experimental $J_{ij}$

In practice, the experimental values of the long-range coupling are not
completely described by a power-law decay, and at large distances decay
exponentially; see section V.2. In spite of this, we show that, based on a
spin-wave analysis, a phase transition still occurs close to the predicted
critical point and that a single soft mode is responsible for the critical
behaviour. To this end, we first note that the experimental Hamiltonian is
given by eq. 1 without the Kac factor. In plotting various quantities however,
we divide by the corresponding Kac factor to compare and contrast different
cases on equal footing.

(a)  Figure S4: The two lowest energy eigenvalues $\Lambda_{0}$ and
$\Lambda_{1}$ versus the (Kac-normalized) transverse field $B/\mathcal{J}$,
corresponding to the experimental $J_{ij}$ for different system sizes;
$\mathcal{J}$ is the Kac factor. The lowest-frequency mode softens,
$\Lambda_{0}\to 0$, approximately at the mean-field predicted critical point
$B_{c}/\mathcal{J}=1$. While the next eigenvalue $\Lambda_{1}$ becomes smaller
for larger and larger system sizes, in contrast with a purely power-law
$J_{ij}$, it remains finite even at the critical point (specifically,
$\Lambda_{1}/{\cal J}\approx 0.5$ for $N=50$ at $B=B_{c}$).

In fig. S4, we plot the first two lowest energy eigenvalues as a function of
$B$ for different system sizes each with the corresponding experimental values
of $J_{ij}$; we have divided $B$ by the Kac factor $\mathcal{J}$ on the
$x$-axis and re-scaled the eigenvalues accordingly. The points where the
lowest-frequency mode softens (i.e., $\Lambda_{0}\to 0)$ agree very well with
the mean-field prediction $B_{c}/\mathcal{J}\sim 1$. However, the gap between
$\Lambda_{1}$ and $\Lambda_{0}$ (in units of the Kac factor) exhibits a strong
dependence on system size, in contrast to the power-law case. This suggests a
departure from the collective nature of the power-law interactions with $p<1$
near the critical point. As $\Lambda_{1}$ becomes smaller, the corresponding
mode becomes populated as well and contributes to the correlations. At the
same time, we note that the latter mode still remains gapped even for $N=50$
where $\Lambda_{1}/{\cal J}\sim 0.5$ at the critical point $B/{\cal J}\approx
1$ while $\Lambda_{0}=0$. Therefore, the soft mode corresponding to
$\Lambda_{0}$ should govern the critical behaviour despite the model’s
departure from infinite- or truly long-range models considered before.

For a quantitative estimate of the contribution due to the first gapped mode
(corresponding to $\Lambda_{1}$), let us first consider a toy model where a
harmonic oscillator in its ground state undergoes a sudden frequency quench
$\Omega_{0}\to\Omega$. This quench excited the harmonic oscillator and creates
a population of $n+1/2=\frac{1}{4}(\Omega/\Omega_{0}+\Omega_{0}/\Omega)$. Now
considering the first gapped mode, we can substitute $\Omega_{0}=2B$ and
$\Omega=\Lambda_{1}$, which gives

$n_{1}\sim\frac{1}{2\Lambda_{1}/{\cal J}}\sim 1\,,$ (56)

where the numerical estimate is computed by setting $B=B_{c}={\cal J}$ and
$\Lambda_{1}/{\cal J}\sim 0.5$ consistent with the data for $N=50$ in fig. S4.
This should be contrasted against the soft mode whose population scales with
system size, according to our scaling analysis, as $n_{0}\sim\sqrt{N}$ at the
critical point. For $N=50$, for example, this means that the soft mode is
populated by $n_{0}\lesssim 10$ in contrast with $n_{1}\sim 1$. While the
higher excited modes may be populated with $n\sim 1$, the soft mode still
governs the critical behaviour at the phase transition.

##  V Experiment vs Theory

In this section, we contrast various theoretical predictions against the
experimental results.

###  V.1 Finite-Size Behavior of the Order Parameter

Mean-field analysis, together with finite-size corrections, yields eq. 28 up
to an overall coefficient. Partial data and the comparison against the
experiment is shown in Fig. 2(b) of the main text. In fig. S5, we provide
additional information for larger system sizes and report the best fit for the
critical point $B_{c}$.

Figure S5: Experimental data for the order parameter against the finite-size
corrected mean-field prediction in eq. 28; $\cal J$ is the Kac normalization.
The experimental data are well described by the theoretical finite-size
correction.

We find that the experimental data are in good agreement with the mean-field
prediction $B_{c}/{\cal J}=1$.

###  V.2 Long-Range Character of $J_{ij}$ Matrix

In fig. S6, we present the numerically calculated average interaction
$J(r)=\frac{1}{N-r}\sum_{i}J_{i,i+r}$ as a function of distance $r$ for
$N=10$, and $50$. We first fit the numerical values to a pure power-law
$J(r)\sim J(1)/r^{p}$ (dashed line in Fig. S6) and find the exponent $p$ to be
around $0.9$ for both system sizes.

Figure S6: Interaction profile $J(r)$ vs $r$. The numerical values of $J(r)$
are shown for system sizes $N=10,\>50$ (see the squares). A power-law fit,
$J(r)/J(1)\sim r^{-p}$, gives $p=0.89$ for both system sizes (see the dashed
lines), and deviates from $J(r)$ at large distances. A better fit is provided
by the product of a power-law with an exponential decay, $J(r)/J(1)\sim
r^{-p}e^{-k(r-1)}$, and gives $p=0.306,\>k=0.231$ for $N=10$, and
$p=0.306,\>k=0.135$ for $N=50$ (see the solid lines). The latter fit captures
the trend of $J(r)$ at larger distances as well.

However, we note that the power-law fit deviates quite significantly from the
calculated values, specifically at large distances and for larger system
sizes. We further fit $J(r)$ against the product of a power-law and an
exponentially decaying function of the form
$J(r)=\frac{J(1)}{r^{p}}e^{-k(r-1)}$ Pagano2020 ; see the solid lines in Fig.
S6. We find the exponents $p=0.306,\>k=0.231$ for $N=10$ and
$p=0.306,\>k=0.135$ for $N=50$ and note that the fit represents the numerical
values closely. The similarity of the power-law exponents for two system sizes
indicates that the experimental interaction matrices are self-similar at short
distances. In section V.3, we show that the spatial correlations after a
single quench shows a slower decay compared to the spin-spin interactions.
This behaviour is consistent with the critical correlations building up in the
system.

###  V.3 Correlation Matrix

In this section, we provide the correlation matrix
$\left\langle\sigma_{i}^{x}\sigma_{j}^{x}\right\rangle$ after a single quench.
Specifically, we consider the average correlations between sites at a fixed
distance $r$ and at the same time as the peak of $\left\langle
S_{x}^{2}\right\rangle$, defined as
${C(r)=\frac{1}{N-r}\sum_{j}\left\langle\sigma_{j}^{x}\sigma^{x}_{j+r}\right\rangle}$.
In fig. S7, we show the correlations in the experiment and contrast them
against the exact simulation of the experimental model when available. Figure
S7(a) shows that the correlation function is in reasonable agreement with the
exact simulation. Sources of the remaining discrepancy are expected to include
imperfect experimental detection and decoherence effects.

Figure S7(b) shows the correlation function for a system of size $N=50$.
Again, we find a much slower decay of the correlation function compared to the
interaction profile $J(r)$, indicating criticality and long-range correlations
in the system. On the other hand, the decay of the correlation function is
relatively pronounced, a feature that should be contrasted against the LMG
model (see section II) which shows no decay due to infinite-range
interactions, and the purely power-law models (see section III) which are
expected to show a relatively slow decay of the correlations with distance.
This observation once again underscores the departure of the experiment from
infinite- or truly long-range models with $p<1$. In light of this, it is quite
remarkable that the critical behaviour observed in the experiment is almost
identical to the latter models (cf. section IV).

Figure S7: Correlations upon a sudden quench to the critical point for (a)
$N=15$ and (b) $N=50$; these correlations are measured/computed at the same
time as the peak of $\left\langle S_{x}^{2}\right\rangle$. (a) Correlations
decay slowly compared to $J(r)$ indicating the build-up of critical
fluctuations. (b) The data for $N=50$ shows consistent results where
correlations decay rather slowly compared to the interaction $J(r)$,
indicating critical correlations. On the other hand, correlations decay
significantly with $r$ in contrast with the infinite- or truly long-range
chain chains, reflecting the distinct character of the experimental system
(see section V.2).

###  V.4 Double Quench Switch Times

Our experimental and simulation data for the critical-to-critical quench does
not collapse perfectly, even with exponents found by optimizing our objective
function. We have found this to be primarily due to deviations in the time at
which the second quench was performed in the experimental system, which were
then used in the simulations.

The second quench should be performed at the same characteristic time for all
system sizes; that is, the same critically scaled time of the first quench,
$\mathcal{J}t/N^{\zeta_{1}}$ where $\zeta_{1}\approx 1/4$. However, in order
to avoid using the expected results in the data analysis, in Fig. 4 of the
main text we restricted ourselves to choosing the switch time for each dataset
(whether experiment or simulation) as the time at which the fluctuations in
the experimental dataset were observed to reach a maximum, which has an
inherent uncertainty due to the limited data sampling and experimental error
sources. As a result, the switch times deviated around this value between
system sizes.

Figure S8: Simulation of the experimental critical-to-critical quench using
corrected switch times. These are to be compared to Fig. 4b,d of the main
text, which instead use the best experimental estimates for these switch
times. The critical exponents are calculated using the dynamics up to the
first minimum after the first maximum, just as in the main text. The collapse
up to and including the first peak is nearly perfect with the critical
exponents that are in much better agreement with the theory prediction.

In fig. S8, we show results from simulating the experimental system with
identical characteristic times. This figure demonstrates that, with the
corrected quench times, the data collapse nearly exactly.

##  VI Optimal Collapse Algorithm

To determine the scaled data shown in Figs 3c,d and 4c,d of the main text, we
construct an objective function to measure the quality of the collapse of
several curves each with some uncertainty. Previous work has used an objective
function similar to the weighted-squared-difference, where all curves are
linearly interpolated together to produce an estimate for the so-called
“master-curve,” from which the error is calculated as the distance to each
curve Houdayer04 . Here we propose an alternative objective function that is
more general and less sensitive to noise in the input data.

We define a curve in this context as a set of data points with associated
uncertainties. Each point is denoted by $\mathbf{p}\equiv[x,y,\Delta x,\Delta
y]\equiv[\mathbf{p}_{x},\mathbf{p}_{y},\Delta\mathbf{p}_{x},\Delta\mathbf{p}_{y}]$,
and a curve with $k$ data points by
$C\equiv\\{\mathbf{p}\\}\in\mathbb{R}^{4\times k}$. We first consider the case
of two curves, $C$ and $C^{\prime}$ with $k$ and $k^{\prime}$ data points with
$k$ and $k^{\prime}$ not necessarily equal. We begin by equally normalizing
the two curves in both $x$ and $y$ such that they are wholly contained within
the region $[0,1]\times[0,1]$. Denote these normalized curves by
$\overline{C}$ and $\overline{C^{\prime}}$ given by

$\displaystyle\begin{aligned} Y_{\mathrm{max}}&=\max_{\mathbf{p}\in C\cup
C^{\prime}}\mathbf{p}_{y},\\\ Y_{\mathrm{min}}&=\min_{\mathbf{p}\in C\cup
C^{\prime}}\mathbf{p}_{y},\end{aligned}$ $\displaystyle\begin{aligned}
X_{\mathrm{max}}&=\max_{\mathbf{p}\in C\cup C^{\prime}}\mathbf{p}_{x},\\\
X_{\mathrm{min}}&=\min_{\mathbf{p}\in C\cup
C^{\prime}}\mathbf{p}_{x},\end{aligned}$ (57)
$\displaystyle\overline{C}=\left\\{\frac{\mathbf{p}-X_{\mathrm{min}}}{X_{\mathrm{max}}-X_{\mathrm{min}}},\frac{\mathbf{p}_{y}-Y_{\mathrm{min}}}{Y_{\mathrm{max}}-Y_{\mathrm{min}}},\right.$
$\displaystyle\left.\frac{\Delta\mathbf{p}_{x}}{X_{\mathrm{max}}-X_{\mathrm{min}}},\frac{\Delta\mathbf{p}_{y}}{Y_{\mathrm{max}}-Y_{\mathrm{min}}}\middle|\mathbf{p}\in
C\right\\},$
$\displaystyle\overline{C^{\prime}}=\left\\{\frac{\mathbf{p^{\prime}}-X_{\mathrm{min}}}{X_{\mathrm{max}}-X_{\mathrm{min}}},\frac{\mathbf{p^{\prime}}_{y}-Y_{\mathrm{min}}}{Y_{\mathrm{max}}-Y_{\mathrm{min}}},\right.$
$\displaystyle\left.\frac{\Delta\mathbf{p^{\prime}}_{x}}{X_{\mathrm{max}}-X_{\mathrm{min}}},\frac{\Delta\mathbf{p^{\prime}}_{y}}{Y_{\mathrm{max}}-Y_{\mathrm{min}}}\middle|\mathbf{p^{\prime}}\in
C^{\prime}\right\\}.$

We define the loss function for a single point in $\overline{C}$ relative to
$\overline{C^{\prime}}$,
$L\left(\mathbf{p}\in\overline{C},\overline{C^{\prime}}\right)$, using the
square of the distance from the point to the nearest line segment in
$\overline{C^{\prime}}$, weighted by the square of the uncertainty in that
distance. If the point in $\overline{C}$ falls outside of the bounds of
$\overline{C^{\prime}}$ then the loss is $0$.

$L\left(\mathbf{p}\in\overline{C},\overline{C^{\prime}}\right)=\begin{cases}0&\text{if
}\mathbf{p}_{x}<\min\limits_{\mathbf{p^{\prime}}\in
C^{\prime}}\mathbf{p^{\prime}}_{x}\text{ or
}\mathbf{p}_{x}>\max\limits_{\mathbf{p^{\prime}}\in
C^{\prime}}\mathbf{p^{\prime}}_{x},\\\
D^{2}\left(\mathbf{p},\overline{C^{\prime}}\right)/\Delta^{2}\left(\mathbf{p},\overline{C^{\prime}}\right)&\text{otherwise},\end{cases}$
(58)

where $D^{2}\left(\mathbf{p},\overline{C^{\prime}}\right)$ is the square of
the distance from $\mathbf{p}$ to the nearest line segment in
$\overline{C^{\prime}}$,

$\displaystyle D^{2}\left(\mathbf{p},\overline{C^{\prime}}\right)$
$\displaystyle\equiv\frac{\left[\left(\mathbf{q^{\prime}}_{x}-\mathbf{r^{\prime}}_{x}\right)\left(\mathbf{r^{\prime}}_{y}-\mathbf{p}_{y}\right)-\left(\mathbf{r^{\prime}}_{x}-\mathbf{p}_{x}\right)\left(\mathbf{q^{\prime}}_{y}-\mathbf{r^{\prime}}_{y}\right)\right]^{2}}{\left(\mathbf{q^{\prime}}_{x}-\mathbf{r^{\prime}}_{x}\right)^{2}+\left(\mathbf{q^{\prime}}_{y}-\mathbf{r^{\prime}}_{y}\right)^{2}},$
(59) $\displaystyle\mathbf{q^{\prime}}$
$\displaystyle\equiv\begin{subarray}{c}\mathrm{argmin}\\\
\mathbf{b^{\prime}}\in\overline{C^{\prime}}\\\
\mathbf{b^{\prime}}_{x}\geq\mathbf{p}_{x}\end{subarray}\left|\mathbf{b^{\prime}}_{x}-\mathbf{p}_{x}\right|,$
$\displaystyle\mathbf{r^{\prime}}$
$\displaystyle\equiv\begin{subarray}{c}\mathrm{argmin}\\\
\mathbf{b^{\prime}}\in\overline{C^{\prime}}\\\
\mathbf{b^{\prime}}_{x}\leq\mathbf{p}_{x}\end{subarray}\left|\mathbf{b^{\prime}}_{x}-\mathbf{p}_{x}\right|,$

and $\Delta^{2}\left(\mathbf{p},\overline{C^{\prime}}\right)$ is the square of
the uncertainty in that distance, found by propagating the uncertainties in
the relevant parameters. Note that points which fall outside the domain of
$C^{\prime}$ are given a loss of $0$, since they have no corresponding line
segment in $C^{\prime}$. A diagram of this calculation is shown in fig. S9.

Figure S9:  Illustration of the calculations involved in the loss function,
eq. 58.

From this formulation it is clear that the loss of a point in $\mathbf{p}\in
C$ that falls exactly on some part of $C^{\prime}$ is $0$, and that if the
uncertainty of $C^{\prime}$ in the region of $\mathbf{p}$ is large then the
loss will be small.

We can then define the cost function between the two curves $C$ and
$C^{\prime}$, $s\left(C,C^{\prime}\right)$ as the sum of the losses of each
point in $\overline{C}$ relative to $\overline{C^{\prime}}$, and vice versa,
normalized by the number of points which were not excluded in eq. 58,

$s\left(C,C^{\prime}\right)\equiv\frac{1}{2}\left[\frac{\sum\limits_{\mathbf{p}\in\overline{C}}L\left(\mathbf{p},\overline{C^{\prime}}\right)}{\left|\overline{C}\right|_{\overline{C^{\prime}}}}+\frac{\sum\limits_{\mathbf{p^{\prime}}\in\overline{C^{\prime}}}L\left(\mathbf{p^{\prime}},\overline{C}\right)}{\left|\overline{C^{\prime}}\right|_{\overline{C}}}\right],$
(60)

where

$\left|A\right|_{B}\equiv\sum\limits_{\begin{subarray}{c}\mathbf{a}\in A\\\
\min\limits_{\mathbf{b}\in
B}\mathbf{b}_{x}\leq\mathbf{a}_{x}\leq\max\limits_{\mathbf{b}\in
B}\mathbf{b}_{x}\end{subarray}}$ (61)

denotes the number of elements in $A$ which fall within the domain of $B$.

Finally, the objective function for determining the collapse of a set of $N$
curves $C_{1},C_{2},\ldots,C_{N}$ is

$S(C_{1},C_{2},\ldots,C_{N})\equiv\frac{1}{N^{2}-N}\sum\limits_{i<j}^{N}s\left(C_{i},C_{j}\right).$
(62)

The normalization of $N^{2}-N$ arises from the number of pairs of non-
identical curves in the set.

We perform convex optimization on this objective function with our data using
Powell’s method Powell64  in order to determine the critical exponents which
result in an optimal collapse. Figure S10 shows the output space of our
objective function, demonstrating that it is well-behaved and exhibits an
identifiable minimum.

(a)

(b)

Figure S10:  Output space of the objective function, eq. 62, around the
calculated minimum for the (a) simulated and (b) experimental type-I quenches.

###  VI.1 Uncertainty Estimates via the Hessian

Equation 62 incorporates the uncertainty in the data directly. To obtain an
uncertainty in the predicted minimum location of the objective function, we
use the Hessian. The uncertainty of any quantity $X$ as a function of the
input parameters $\\{x_{i}\\}$ of an objective function $f$ can be estimated
by Pumplin01

$\Delta\left(X\right)^{2}=\left(\tilde{f}-\hat{f}\right)\sum\limits_{i,j}H_{ij}^{-1}\frac{\partial
X}{\partial x_{i}}\frac{\partial X}{\partial x_{j}}.$ (63)

Where $H^{-1}$ is the inverse of the Hessian matrix of $f$ evaluated at the
location of the predicted minimum, $\tilde{f}$ is the value of $f$ evaluated
at the same location, and $\hat{f}$ is the true global minimum value of the
objective function.

In our application $\\{x_{i}\\}=\\{\alpha,\zeta\\}$, and $X=\alpha$ or $\zeta$
as well. Moreover, we have designed our objective function such that the
global minimum occurs when the curves perfectly collapse and so $\hat{f}=0$.
Therefore the uncertainty in the predicted exponents $\tilde{\alpha}$ and
$\tilde{\zeta}$ is simply given by

$\displaystyle\Delta\tilde{\alpha}$
$\displaystyle=\sqrt{S\left(\tilde{\alpha},\tilde{\zeta}\right)H_{\alpha\alpha}^{-1}},$
(64) $\displaystyle\Delta\tilde{\zeta}$
$\displaystyle=\sqrt{S\left(\tilde{\alpha},\tilde{\zeta}\right)H_{\zeta\zeta}^{-1}}.$

## References

  * (1) Islam R., Senko C., Campbell W. C., Korenblit S., Smith J., Lee A., Edwards E. E., Wang C.-C. J., Freericks J. K., and Monroe C.  Emergence and frustration of magnetism with variable-range interactions in a quantum simulator.  Science, 340(6132):583–587, 2013.
  * (2) Jurcevic P., Shen H., Hauke P., Maier C., Brydges T., Hempel C., Lanyon B. P., Heyl M., Blatt R., and Roos C. F.  Direct observation of dynamical quantum phase transitions in an interacting many-body system.  Phys. Rev. Lett., 119:080501, Aug 2017.
  * (3) Tan W. L., Becker P., Liu F., Pagano G., Collins K., De A., Feng L., Kaplan H., Kyprianidis A., Lundgren R., et al.  Domain-wall confinement and dynamics in a quantum simulator.  Nature Physics, 17(6):742–747, 2021.
  * (4) Suzuki S., Inoue J.-i., and Chakrabarti B. K.  Quantum Ising phases and transitions in transverse Ising models, volume 862.  Springer, 2012.
  * (5) Titum P. and Maghrebi M. F.  Nonequilibrium criticality in quench dynamics of long-range spin models.  Phys. Rev. Lett., 125:040602, Jul 2020.
  * (6) Defenu N.  Metastability and discrete spectrum of long-range systems.  Proceedings of the National Academy of Sciences, 118(30):e2101785118, 2021.
  * (7) Hauke P. and Tagliacozzo L.  Spread of correlations in long-range interacting quantum systems.  Phys. Rev. Lett., 111:207202, Nov 2013.
  * (8) Colpa J. H. P.  Diagonalization of the quadratic boson hamiltonian.  Physica A: Statistical Mechanics and its Applications, 93(3):327–353, September 1978.
  * (9) Pagano G., Bapat A., Becker P., Collins K. S., De A., Hess P. W., Kaplan H. B., Kyprianidis A., Tan W. L., Baldwin C., Brady L. T., Deshpande A., Liu F., Jordan S., Gorshkov A. V., and Monroe C.  Quantum approximate optimization of the long-range Ising model with a trapped-ion quantum simulator.  PNAS, 117(41), (2020).
  * (10) Houdayer J. and Hartmann A. K.  Low-temperature behavior of two-dimensional gaussian ising spin glasses.  Phys. Rev. B, 70:014418, Jul 2004.
  * (11) Powell M. J. D.  An efficient method for finding the minimum of a function of several variables without calculating derivatives.  The Computer Journal, 7(2):155–162, 01 1964.
  * (12) Pumplin J., Stump D. R., and Tung W. K.  Multivariate fitting and the error matrix in global analysis of data.  Phys. Rev. D, 65:014011, Dec 2001.