# Large Deviations Theory Beyond the Kibble-Zurek Mechanism

Federico Balducci Department of Physics and Materials Science, University of
Luxembourg, L-1511 Luxembourg, G. D. Luxembourg  Mathieu Beau Department of
Physics, University of Massachusetts, Boston, MA 02125, USA  Jing Yang
Department of Physics and Materials Science, University of Luxembourg, L-1511
Luxembourg, G. D. Luxembourg  Nordita, KTH Royal Institute of Technology and
Stockholm University, Hannes Alfvéns väg 12, SE-106 91 Stockholm, Sweden
Andrea Gambassi SISSA - International School for Advanced Studies, via Bonomea
265, 34136 Trieste, Italy  INFN, Sezione di Trieste, Trieste, Italy  Adolfo
del Campo<EMAIL_ADDRESS>Department of Physics and Materials Science,
University of Luxembourg, L-1511 Luxembourg, G. D. Luxembourg  Donostia
International Physics Center, E-20018 San Sebastián, Spain

###### Abstract

The Kibble-Zurek mechanism (KZM) predicts that the average number of
topological defects generated upon crossing a quantum phase transition obeys a
universal scaling law with the quench time. Fluctuations in the defect number
near equilibrium are approximately of Gaussian form, in agreement with the
central limit theorem. Using large deviations theory, we characterize the
universality of fluctuations beyond the KZM and report the exact form of the
rate function in the transverse-field quantum Ising model. In addition, we
characterize the scaling of large deviations in an arbitrary continuous phase
transition, building on recent evidence establishing the universality of the
defect number distribution.

The Kibble-Zurek mechanism (KZM) is an important paradigm in nonequilibrium
statistical physics, describing the dynamics across a continuous phase
transition Polkovnikov  _et al._ (2011); Dziarmaga (2010); del Campo and Zurek
. The divergence of the equilibrium relaxation time in the neighborhood of the
critical point makes the critical dynamics necessarily nonadiabatic for large
systems and leads to the spontaneous formation of topological defects.
Consider a phase transition from a high symmetry phase to a broken symmetry
phase, induced by varying a control parameter $g$ across its critical value
$g_{c}$ in a finite quench time $\tau_{Q}$. The central prediction of the KZM
is that the average defect density, generated during the phase transition,
displays a universal power-law dependence as a function of the quench time.
The KZM thus makes a quantitative prediction on the breakdown of adiabatic
dynamics across a phase transition and holds both in the classical and quantum
regimes Polkovnikov (2005); Zurek  _et al._ (2005); Damski (2005); Dziarmaga
(2005); Polkovnikov  _et al._ (2011).

Kibble’s pioneering work was motivated by cosmological considerations
regarding structure formation in the early universe Kibble (1976). The
prospect of exploring analogous phenomena in condensed matter systems was soon
realized Kibble (1980); Zurek (1985, 1996) and pursued experimentally del
Campo and Zurek ; Deutschländer  _et al._ (2015); Maegochi  _et al._ (2022);
Du  _et al._ (2023). The advance of quantum technologies has led to new tests
of the KZM using quantum simulators in a variety of platforms, including
ultracold gases Weiler  _et al._ (2008); Lamporesi  _et al._ (2013); Navon
_et al._ (2015); Anquez  _et al._ (2016); Ko  _et al._ (2019); Yi  _et al._
(2020); Qiu  _et al._ (2020), trapped ions Ulm  _et al._ (2013); Pyka  _et
al._ (2013); Cui  _et al._ (2016, 2020); Li  _et al._ (2023), and Rydberg
gases Keesling  _et al._ (2019); Ebadi  _et al._ (2021). Recently, the KZM has
been studied with quantum computing devices, such as quantum annealers Gardas
_et al._ (2018); Weinberg  _et al._ (2020); Bando  _et al._ (2020); King  _et
al._ (2022). The accumulated body of literature broadly supports the validity
of KZM in a wide variety of systems.

Experiments probing critical dynamics generally involve an ensemble of single
experimental runs or individual realizations in which measurements are
performed. As a result, they can access information beyond the average defect
density and characterize the ensemble statistics. It is thus natural to ask
whether there are universal signatures in the statistical properties of
spontaneously generated topological defects Cincio  _et al._ (2007); del Campo
(2018); del Campo  _et al._ (2022). The full counting statistics of defects
appears to be universal in classical and quantum systems. Specifically, in
classical continuous phase transitions, it has been found that the defect
number distribution is binomial with an average density in agreement with the
KZM Gómez-Ruiz  _et al._ (2020); Mayo  _et al._ (2021); del Campo  _et al._
(2021); Gómez-Ruiz  _et al._ (2022). Exact solutions in quantum integrable
systems have shown that the kink number distribution is Poisson-binomial del
Campo (2018); Cui  _et al._ (2020), a feature that can hold even when the
system is coupled to an environment Bando  _et al._ (2020); King  _et al._
(2022). These predictions build on the conventional KZM but lie outside its
scope, requiring additional assumptions. We shall thus refer to them as
beyond-KZM physics.

The average number of defects is an extensive quantity. By contrast, the
defect density is intensive, and its fluctuations near equilibrium are
approximately Gaussian, in agreement with the central limit theorem. Large
deviations theory (LDT) addresses the probability of nontypical events in
which an intensive quantity deviates from its average value. The probability
of such large deviations decays exponentially upon increasing the system size,
with a rate controlled by the so-called rate function Ellis (2006); Touchette
(2009). LDT provides a building block of statistical mechanics in and out of
equilibrium. As such, it is a natural framework to explore beyond-KZM physics.
To date, LDT has been used to describe the dynamics of many-body quantum
systems in the limit of sudden quenches when $\tau_{Q}\rightarrow 0$, e.g., to
characterize the work statistics of a given process Gambassi and Silva (2012);
Goold  _et al._ (2018); Perfetto  _et al._ (2019).

In this Letter, we establish the universality of large deviations beyond the
KZM, after crossing a quantum phase transition in a finite time. Specifically,
we report the exact rate function of the driven transverse-field quantum Ising
model (TFQIM), characterizing the statistics of large fluctuations away from
the mean kink density predicted by the conventional KZM. We further generalize
these results to characterize the universality of large deviations in an
arbitrary continuous phase transition leading to point-like defects.

The transverse-field quantum Ising model. The TFQIM has been instrumental in
generalizing the KZM from the classical to the quantum domain Polkovnikov
(2005); Zurek  _et al._ (2005); Damski (2005); Dziarmaga (2005); Cui  _et al._
(2016), and assessing the universality of beyond-KZM physics, both in theory
del Campo (2018) and experiments Cui  _et al._ (2020); Bando  _et al._ (2020);
King  _et al._ (2022). Its Hamiltonian is given by

$H[g(t)]=-J\sum_{l=1}^{N}\left[g(t)\sigma_{l}^{x}+\sigma_{l}^{z}\sigma_{l+1}^{z}\right],$
(1)

where $J>0$ favors ferromagnetic alignment and $g(t)$ plays the role of an
effective magnetic field. In the fermionic representation, the Ising chain
Hamiltonian becomes Sup

$H[g(t)]=2J\sum_{k>0}\psi_{k}^{\dagger}\left[\tau^{z}(g(t)-\cos
k)+\tau^{y}\sin k\right]\psi_{k},$ (2)

in terms of the fermionic operators
$\psi_{k}^{\dagger}\equiv(\tilde{c}_{k}^{\dagger},\tilde{c}_{-k}^{\phantom{\dagger}})$
in momentum space. Here, $\tau^{x,y,z}$ are Pauli matrices. We choose to work
with periodic boundary conditions so that the momentum is a good quantum
number and takes the values $k=(2n+1)\pi/N$ with $n=-N/2,\dots,N/2-1$, as
discussed, e.g., in Refs. Damski and Rams (2013); Sup . Momentum conservation
restricts the formation of defects to kink-antikink pairs. Choosing the total
number of spins $N$ to be even proves convenient since the number of kink
pairs is then restricted to outcomes in the set $\\{0,1,2,\dots,N/2\\}$. Given
Eq. (2), the dynamics of the TFQIM can be reduced to that of an ensemble of
non-interacting two-level systems Dziarmaga (2005).

Consider a quench, in a finite time $\tau_{Q}$, from the paramagnetic to the
ferromagnetic phase

$g(t)=g_{c}\left(1-\frac{t}{\tau_{Q}}\right),$ (3)

where $g(0)=g_{c}=1$ is the critical value of $g$, and we let $t$ run from
$-3\tau_{Q}$ to $\tau_{Q}$. We will refer to $\tau_{Q}$ as the quench time. We
choose $g(\tau_{Q})=0$ for simplicity since the final Hamiltonian contains
only the ferromagnetic term and commutes with the kink-pair number operator
$K_{N}\equiv\frac{1}{4}\sum_{l=1}^{N}\left(1-\sigma^{z}_{l}\sigma^{z}_{l+1}\right)$.
This observable counts the number of kink-antikink pairs in a given quantum
state and is extensive in the system size $N$ Dziarmaga (2005). The study of
its eigenvalue statistics provided the basis of previous studies exploring
universality beyond the KZM del Campo (2018); Cui  _et al._ (2020); Bando  _et
al._ (2020); Mayo  _et al._ (2021); King  _et al._ (2022). We define an
intensive kink-pair density operator

$\hat{\rho}_{N}\equiv\frac{K_{N}}{N}=\frac{1}{4N}\sum_{l=1}^{N}\left(1-\sigma^{z}_{l}\sigma^{z}_{l+1}\right).$
(4)

The density of kink pairs, upon completion of the quench in Eq. (3), is given
by the expectation value $\rho_{\rm KZM}=\langle\hat{\rho}_{N}\rangle$ at the
final time $\tau_{Q}$. It exhibits a power-law scaling in the slow driving
limit, i.e., to leading order in a $1/\tau_{Q}$ expansion Polkovnikov (2005);
Zurek  _et al._ (2005); Damski (2005); Dziarmaga (2005); Polkovnikov  _et al._
(2011)

$\rho_{\rm
KZM}=\langle\hat{\rho}_{N}\rangle=\frac{1}{4\pi}\sqrt{\frac{\hbar}{2J\tau_{Q}}},$
(5)

in agreement with the celebrated, universal KZM power-law scaling $\rho_{\rm
KZM}\propto\tau_{Q}^{-\frac{\nu}{1+z\nu}}$ for the critical exponents
$\nu=z=1$ of the TFQIM del Campo and Zurek .

In any quantum state other than an eigenstate of $H(g=0)$, the density
operator $\hat{\rho}_{N}$ will exhibit fluctuations of either classical or
quantum nature. The probability distribution function $P(\rho_{N})$,
characterizing the eigenvalue statistics of the kink-pair density operator,
reads

$P(\rho_{N})=\left\langle\delta(\hat{\rho}_{N}-\rho_{N})\right\rangle,$ (6)

where $\rho_{N}$ is the random variable associated with the kink-pair-number
operator $\hat{\rho}_{N}$. We aim at uncovering via LDT the universality of
large fluctuations of $P(\rho_{N})$ away from the mean, which the conventional
KZM predicts.

Large deviations theory beyond the KZM in the TFQIM. The central object in LTD
is the scaled cumulant generating function, associated with a random variable
$\rho_{N}$, depending on a large parameter $N$,

$\lambda(\theta)=\lim_{N\rightarrow\infty}\frac{1}{N}\ln\left\langle
e^{N\theta\hat{\rho}_{N}}\right\rangle.$ (7)

The Gärtner-Ellis theorem states that when $\lambda(\theta)$ exists for all
$\theta\in\mathbb{R}$, then the random variable $\rho_{N}$ satisfies the large
deviations principle Ellis (2006); Touchette (2009)

$P\big{(}\rho_{N}\in[\rho,\rho+d\rho]\big{)}\approx e^{-NI(\rho)}d\rho,$ (8)

with the rate function $I(\rho)$ given by the Legendre-Fenchel transform

$I(\rho)=\sup_{\theta\in\mathbb{R}}\,\big{[}\theta\rho-\lambda(\theta)\big{]}.$
(9)

Deviations from the mean value are thus exponentially suppressed by the rate
function $I(\rho)$ weighted with the system size, and the random variable
concentrates around the mean in the thermodynamic limit.

Let us consider the application of the Gärtner-Ellis theorem to the
distribution of kink pairs generated across a quantum phase transition in the
TFQIM. In this case, the defect density is a non-negative quantity. As a
result, $I(\rho)$ is divergent for $\rho<0$, and we focus on the domain
$\rho\geq 0$. We note that in Fourier space, the operator associated with the
density of kink pairs at the end of the quench is

$\hat{\rho}_{N}=\frac{1}{N}\sum_{k>0}\gamma_{k}^{\dagger}(\tau_{Q})\gamma_{k}^{\phantom{\dagger}}(\tau_{Q}),$
(10)

where $\gamma_{k}^{\phantom{\dagger}}(\tau_{Q})$ and
$\gamma_{k}^{\dagger}(\tau_{Q})$ are the fermionic Bogoliubov operators at the
end of the quench, and the sum is restricted to $k>0$ since the number of kink
pairs equals the number of right-moving kinks. Further, for free fermions
(with periodic boundary conditions), the time-dependent density matrix
$\varrho(t)$ retains the tensor product structure during unitary time
evolution, i.e., $\varrho(t)=\bigotimes_{k}\varrho_{k}(t)$. As a result, the
moment-generating function admits the explicit form

$\displaystyle\left\langle e^{N\theta\hat{\rho}_{N}}\right\rangle$
$\displaystyle=\prod_{k>0}\operatorname{Tr}\left[\varrho_{k}(\tau_{Q})e^{\theta\gamma_{k}^{\dagger}(\tau_{Q})\gamma_{k}^{\phantom{\dagger}}(\tau_{Q})}\right]$
$\displaystyle=\prod_{k>0}\left[1+\left(e^{\theta}-1\right)p_{k}\right],$ (11)

where
$p_{k}=\langle\gamma_{k}^{\dagger}(\tau_{Q})\gamma_{k}^{\phantom{\dagger}}(\tau_{Q})\rangle\in[0,1]$
represents the probability that the mode $k$ is excited at the end of the
protocol. This is the moment-generating function of a Poisson binomial
distribution associated with the sum of $N/2$ independent random Bernoulli
variables, each of which has probability $p_{k}$ for the occupation number to
be $1$, corresponding to the formation of a kink-antikink pair, and
probability $(1-p_{k})$ for the occupation number to be $0$, corresponding to
no defect formation del Campo (2018). In addition, the value of $p_{k}$ can be
estimated according to the Landau-Zener approximation Dziarmaga (2005),
$p_{k}=\langle\gamma_{k}^{\dagger}(\tau_{Q})\gamma_{k}^{\phantom{\dagger}}(\tau_{Q})\rangle\approx\exp(-2\pi
J\tau_{Q}k^{2}/\hbar)$ near $k=0$, dictating an exponential decay with
increasing quench time and a Gaussian decay as a function of the wavenumber.
This behavior dictates the KZM scaling in a quantum phase transition Dziarmaga
(2005); Damski (2005); Damski and Zurek (2006). The explicit computation of
the scaled cumulant generating function, according to Eqs. (7) and (11), in
the limit $N\to\infty$, yields

$\lambda(\theta)=\int_{0}^{\pi}\frac{dk}{2\pi}\,\ln\left[1+\left(e^{\theta}-1\right)p_{k}\right],$
(12)

which is a convergent integral.

For slow quenches, using a power-series expansion in $1/\tau_{Q}$ to leading
order, or equivalently extending the upper limit of the integral in Eq. (12)
to infinity, one finds

$\lambda(\theta)=-\rho_{\rm KZM}\mathrm{Li}_{3/2}\left(1-e^{\theta}\right),$
(13)

in terms of the polylogarithm function ${\rm
Li}_{q}(x)=\sum_{s=1}^{\infty}x^{s}/s^{q}$. The exact expression in the slow-
quench limit, Eq. (13), shows that $\lambda(\theta)$ is differentiable for all
values of $\theta$, ensuring the applicability of the Gärtner-Ellis theorem.
We verify that for $\theta=0$, $\lambda(0)=0$, consistently with its
definition. Further, for $\theta<0$, $\lambda(\theta)$ quickly approaches the
constant value $\lambda(-\infty)=-\rho_{\rm KZM}\zeta(3/2)$. Indeed,
$\lambda(\theta)$ is approximately constant for $\theta<0$ and is a monotonic
function of $\theta$.

We define a dimensionless density of defects $\bar{\rho}\equiv\rho/\rho_{\rm
KZM}$ in terms of which

$I(\rho)=\rho_{\mathrm{KZM}}\sup_{\theta\in\mathbb{R}}\left[\theta\bar{\rho}+{\rm
Li}_{3/2}\left(1-e^{\theta}\right)\right].$ (14)

As a result, the rate function (14) is universal in the sense that
$\bar{I}(\bar{\rho})=I(\rho)/\rho_{\mathrm{KZM}}$ varies only with
$\bar{\rho}$ and is independent of the quench time $\tau_{Q}$. This is the
central result of our work, which we elaborate and generalize in what follows.
Taking the derivative with respect to $\theta$, one finds at the supremum
$\theta^{*}$

$\bar{\rho}=-\frac{e^{\theta^{*}}}{e^{\theta^{*}}-1}\mathrm{Li}_{1/2}\left(1-e^{\theta^{*}}\right).$
(15)

The function $\theta^{*}(\bar{\rho})$ and the rate function scaled by the KZM
density $\bar{I}(\bar{\rho})$ are found numerically. The rate function is
shown in Fig. 1. As the decay of the probability density function $P(\rho)$ is
dictated by the rate function according to Eq. (8), the minimum of
$\bar{I}(\bar{\rho})$ at $\bar{\rho}=1$ is associated with the most likely
value of $\hat{\rho}_{N}$, which equals the mean value $\rho_{\rm KZM}$
predicted by the KZM. Thus, LDT guarantees that the KZM prediction holds with
maximum probability.

Figure 1: Comparison of the scale rate function
$\bar{I}(\bar{\rho})=I(\rho)/\rho_{\mathrm{KZM}}$ derived analytically with
the numerically exact computation for a finite $\tau_{Q}$ and $N$. A TFQI
chain is initialized in its ground state $g(-3\tau_{Q})=4g_{c}$. The equations
of motion are integrated numerically up to time $\tau_{Q}$, corresponding to
$g(\tau_{Q})=0$. From the final wavefunction, the cumulant generating function
$\lambda(\theta)$ is computed using Eq. (7) for finite $N$, from which the
scaled rate function $\bar{I}$ is found with a Legendre-Fenchel transform. As
the quench time increases, the agreement between numerics and the analytical
solution based on the LZ approximation improves visibly, while the agreement
with the central limit theorem (CLT) prediction does not. The value at the
origin is $\bar{I}(0)=\zeta(3/2)$ follows from Eq. (14), while the minimum
$\bar{I}=0$ is attained at $\bar{\rho}=1$ (diamond). Finite-size analysis
reveals the convergence of the numerically-evaluated rate function $\bar{I}$
to the thermodynamic limit for $N=1000$, which is used in this figure.

Concentration inequalities. Let us tackle the problem of large deviations from
a complementary angle using concentration inequalities Vershynin (2018). To
bound the deviations, we make use of the Chernoff bound, which reads
$P(\rho_{N}>\rho)\leq\langle e^{\theta\hat{\rho}_{N}}\rangle e^{-\theta\rho}$,
for all $\theta>0$. The characteristic function can be written as

$\displaystyle\left\langle e^{\theta\hat{\rho}_{N}}\right\rangle$
$\displaystyle=\exp\left\\{\frac{N}{2\pi}\int_{0}^{\pi}dk\,\ln\left[1+(e^{\theta}-1)p_{k}\right]\right\\}$
(16) $\displaystyle\approx\exp\left[-N\rho_{\rm KZM}{\rm
Li}_{3/2}\left(1-e^{\theta}\right)\right].$

We thus find from the Chernoff bound that

$P(\rho_{N}>\rho)\leq\exp\left\\{-\rho_{\rm KZM}[\theta\bar{\rho}+{\rm
Li}_{3/2}(1-e^{\theta})]\right\\}$ (17)

for all $\theta>0$. Tightening the above inequality by taking the supremum of
the exponent, the right tail of the distribution is bounded as

$P(\rho_{N}>\rho)\leq\exp[-NI(\rho)],$ (18)

with $I(\rho)$ given by Eq. (14). Likewise, the left tail is bound by the same
term, $P(\rho_{N}<\rho)\leq\exp[-NI(\rho)]$. The logarithm of the two-sided
Chernoff bound is the rate function. The above results establish the nature of
large deviations of kink-antikink pairs formed across the quantum phase
transition between the paramagnetic and the ferromagnetic phase of the TFQIM.
Such results are generalizable to the family of quasi-free fermion models in
which the density of defects is given by the density of quasiparticles. In
addition, the universal form of the scaled cumulant generating function and
the rate function in the slow quench limit also hold when fast decaying long-
range interactions are considered Sup , further confirming their universality
under fast-decaying long-range deformations. We next turn our attention to an
arbitrary continuous phase transition described by the KZM.

LDT beyond KZM: General scenario. Consider a scenario of spontaneous symmetry
breaking leading to the generation of point-like defects in $d$ spatial
dimensions. The KZM exploits the equilibrium scaling relations for the
correlation length $\xi$ and the relaxation time $\tau$, i.e.,

$\xi=\frac{\xi_{0}}{|\varepsilon|^{\nu}},\qquad\tau=\frac{\tau_{0}}{|\varepsilon|^{z\nu}},$
(19)

where $\nu$ and $z$ are critical exponents and $\xi_{0}$ and $\tau_{0}$ are
microscopic constants. The dimensionless variable
$\varepsilon=(g_{c}-g)/g_{c}$ quantifies the distance to the critical point
$g_{c}$, and vanishes at the phase transition. Linearizing the driving
protocol in the neighbourhood of $g_{c}$ as $g(t)=g_{c}(1-t/\tau_{Q})$, one
identifies the quench time $\tau_{Q}$. The KZM sets the nonequilibrium
correlation length to be
$\hat{\xi}=\xi_{0}(\tau_{Q}/\tau_{0})^{\frac{\nu}{1+z\nu}}$. During the phase
transition, the system is partitioned into protodomains of average volume
$\hat{\xi}^{d}$. A defect may be generated with an empirical probability $p$
at the merging point between adjacent domains. For point-like defects, the
number of events for defect formation is estimated as
$\mathcal{N}=V_{d}/(f\hat{\xi}^{d})$, where $V_{d}$ is the volume in $d$
spatial dimensions and $f$ a fudge factor of order one. As a result, it scales
as
$\mathcal{N}=[V_{d}/(f\xi_{0}^{d})]\left(\tau_{0}/\tau_{Q}\right)^{\frac{d\nu}{1+z\nu}}$.
Assume defect formation events at different locations to be described by
independent and identically distributed discrete random variables $X_{i}$ with
$i=1,\dots,\mathcal{N}$ Gómez-Ruiz  _et al._ (2020); Mayo  _et al._ (2021);
del Campo  _et al._ (2021); Gómez-Ruiz  _et al._ (2022), where the outcome
$X_{i}=+1$ corresponds to the formation of a topological defect, and $X_{i}=0$
to no defect formation. The defect number distribution takes the binomial form
$P(n)=\binom{\mathcal{N}}{n}p^{n}(1-p)^{\mathcal{N}-n}$. Numerical studies
support this prediction in $d=1,2$ for varying $\tau_{Q}$ Gómez-Ruiz  _et al._
(2020); Mayo  _et al._ (2021); del Campo  _et al._ (2021); Gómez-Ruiz  _et
al._ (2022). Accordingly, the average number of topological defects is given
by $\rho_{\rm KZM}=p\mathcal{N}/V_{d}\propto\tau_{Q}^{-\frac{d\nu}{1+z\nu}}$.
The defect density, an intensive random variable, is given by
$\rho_{\mathcal{N}}=\sum_{i=1}^{\mathcal{N}}X_{i}/V_{d}$. We are interested in
estimating the probability distribution of
$S_{\mathcal{N}}=\sum_{i=1}^{\mathcal{N}}X_{i}$ when $\mathcal{N}$ is large.
Using the Stirling approximation, one finds

$P(S_{\mathcal{N}}=r\mathcal{N})=\frac{1}{\sqrt{2\pi
r(1-r)\mathcal{N}}}e^{-V_{d}\rho_{\rm KZM}D_{\rm KL}(r\|p)},$ (20)

where $D_{\rm KL}(r\|p)=r\ln(r/p)+(1-r)\ln[(1-r)/(1-p)]$ is the Kullback-
Liebler distance between two Bernoulli distributions with success
probabilities $r$ and $p$. It satisfies $D_{\rm KL}(r\|p)\geq 0$, with the
equality holding when $r$ equals $p$, which is the most probable value.
Neglecting the prefactor, we thus find that for large $\mathcal{N}$,
fluctuations of the defect number away from the mean are suppressed
exponentially with increasing $\mathcal{N}$, i.e.,
$P(S_{\mathcal{N}}=r\mathcal{N})\approx\exp[-V_{d}\rho_{\rm KZM}D_{\rm
KL}(r\|p)]$. In this sense, the defect number distribution concentrates at the
KZM prediction in the thermodynamic limit when $V_{d}$ and $\mathcal{N}$
diverge. Indeed, in the spirit of LDT, we identify the rate function

$\displaystyle I(r)=\rho_{\rm KZM}D_{\rm KL}(r\|p),$ (21)

generalizing the findings for the TFQIM to arbitrary continuous phase
transitions. Similarly, $I(r)$ dictates the universal suppression of
deviations away from the KZM prediction. For example, the right tail of the
distribution is bounded as

$\displaystyle P(S_{\mathcal{N}}\geq r\mathcal{N})$ $\displaystyle\leq$
$\displaystyle\exp\left[-V_{d}\,I(r)\right]$ $\displaystyle=$
$\displaystyle\exp\left[-\frac{pV_{d}}{f\xi_{0}^{d}}\left(\frac{\tau_{0}}{\tau_{Q}}\right)^{\frac{d\nu}{1+z\nu}}\\!\\!D_{\rm
KL}(r\|p)\right].$

These results are fully consistent with LDT. Using the dimensionless density
of defect $\bar{\rho}\equiv\rho_{\mathcal{N}}/\rho_{\rm
KZM}=\frac{1}{p\mathcal{N}}\sum_{i=1}^{\mathcal{N}}X_{i}$, according to
Sanov’s theorem in LDT Touchette (2009)
$P(p\bar{\rho}=r)=e^{-\mathcal{N}D_{\rm KL}(r\|p)}$. Thus,
$P(\bar{\rho})=\exp[-V_{d}I(\bar{\rho})]$ where $I(\bar{\rho})\equiv\rho_{\rm
KZM}D_{\rm KL}(\bar{\rho}p\|p)/p$. The tails of the distribution read then
$P(\bar{\rho}\geq\sigma)\leq\exp[-V_{d}\,I(\sigma)]$ as in Eq. (18).

Discussion. The rate function governs the nature of deviations away from the
mean according to LDT. Knowledge of the rate function relies on the full
counting statistics of topological defects beyond KZM, as the rate function is
the Legendre-Frenchel transform of the scaled cumulant generating function.
Using the exact solution of the critical dynamics in the TFQIM as a test bed,
we have explored the nature of large deviations in the number of topological
defects generated across a quantum phase transition driven at the finite rate
and show that the rate function is proportional to the KZM density of kink
pairs. The rate function exhibits a universal power-law scaling with the
quench time in which the phase transition is crossed. We have further
generalized these findings to account for the dynamics of arbitrary continuous
phase transitions described by KZM. Assuming the independence of defect
formation events during the phase transition, we have shown that the
probability density function of the defect number distribution decays
exponentially with the number of defects predicted by KZM. Our results are of
broad interest in nonequilibrium quantum and classical statistical mechanics,
connecting LDT with the breakdown of adiabatic dynamics, and should find broad
applications in quantum simulation, quantum annealing, ultracold atom physics,
and the study of critical phenomena.

Acknowledgements. AdC thanks SISSA for its hospitality during the early stages
of this work. We thank Federico Roccati for their feedback on the manuscript.
This project was funded within the QuantERA II Programme that has received
funding from the European Union’s Horizon 2020 research and innovation
programme under Grant Agreement No. 16434093. AG acknowledges financial
support from the PNRR MUR project PE0000023-NQST. For open access and in
fulfillment of the obligations arising from the grant agreement, the authors
have applied a Creative Commons Attribution 4.0 International (CC BY 4.0)
license to any Author Accepted Manuscript version arising from this
submission.

## References

  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Dziarmaga (2010) J. Dziarmaga, Adv. Phys. 59, 1063 (2010).
  * (3) A. del Campo and W. H. Zurek, Intl. J. Mod. Phys. A .
  * Polkovnikov (2005) A. Polkovnikov, Phys. Rev. B 72, 161201 (2005).
  * Zurek  _et al._ (2005) W. H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. 95, 105701 (2005).
  * Damski (2005) B. Damski, Phys. Rev. Lett. 95, 035701 (2005).
  * Dziarmaga (2005) J. Dziarmaga, Phys. Rev. Lett. 95, 245701 (2005).
  * Kibble (1976) T. W. B. Kibble, J. Phys. A: Math. Gen. 9, 1387 (1976).
  * Kibble (1980) T. W. B. Kibble, Phys. Rep. 67, 183 (1980).
  * Zurek (1985) W. H. Zurek, Nature 317, 505 (1985).
  * Zurek (1996) W. H. Zurek, Phys. Rep. 276, 177 (1996).
  * Deutschländer  _et al._ (2015) S. Deutschländer, P. Dillmann, G. Maret, and P. Keim, Proc. the Natl. Acad. Sci. 112, 6925 (2015).
  * Maegochi  _et al._ (2022) S. Maegochi, K. Ienaga, and S. Okuma, Phys. Rev. Lett. 129, 227001 (2022).
  * Du  _et al._ (2023) K. Du, X. Fang, C. Won, C. De, F.-T. Huang, W. Xu, H. You, F. J. Gómez-Ruiz, A. del Campo, and S.-W. Cheong, Nat. Phys. (2023), 10.1038/s41567-023-02112-5.
  * Weiler  _et al._ (2008) C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis, and B. P. Anderson, Nature 455, 948 (2008).
  * Lamporesi  _et al._ (2013) G. Lamporesi, S. Donadello, S. Serafini, F. Dalfovo, and G. Ferrari, Nat. Phys. 9, 656 (2013).
  * Navon  _et al._ (2015) N. Navon, A. L. Gaunt, R. P. Smith, and Z. Hadzibabic, Science 347, 167 (2015).
  * Anquez  _et al._ (2016) M. Anquez, B. A. Robbins, H. M. Bharath, M. Boguslawski, T. M. Hoang, and M. S. Chapman, Phys. Rev. Lett. 116, 155301 (2016).
  * Ko  _et al._ (2019) B. Ko, J. W. Park, and Y. Shin, Nat. Phys. 15, 1227 (2019).
  * Yi  _et al._ (2020) C.-R. Yi, S. Liu, R.-H. Jiao, J.-Y. Zhang, Y.-S. Zhang, and S. Chen, Phys. Rev. Lett. 125, 260603 (2020).
  * Qiu  _et al._ (2020) L.-Y. Qiu, H.-Y. Liang, Y.-B. Yang, H.-X. Yang, T. Tian, Y. Xu, and L.-M. Duan, Sci. Adv. 6, eaba7292 (2020).
  * Ulm  _et al._ (2013) S. Ulm, J. Roßnagel, G. Jacob, C. Degünther, S. T. Dawkins, U. G. Poschinger, R. Nigmatullin, A. Retzker, M. B. Plenio, F. Schmidt-Kaler, and K. Singer, Nat. Commun. 4, 2290 (2013).
  * Pyka  _et al._ (2013) K. Pyka, J. Keller, H. L. Partner, R. Nigmatullin, T. Burgermeister, D. M. Meier, K. Kuhlmann, A. Retzker, M. B. Plenio, W. H. Zurek, A. del Campo, and T. E. Mehlstäubler, Nat. Commun. 4, 2291 (2013).
  * Cui  _et al._ (2016) J.-M. Cui, Y.-F. Huang, Z. Wang, D.-Y. Cao, J. Wang, W.-M. Lv, L. Luo, A. del Campo, Y.-J. Han, C.-F. Li, and G.-C. Guo, Sci. Rep. 6, 33381 (2016).
  * Cui  _et al._ (2020) J.-M. Cui, F. J. Gómez-Ruiz, Y.-F. Huang, C.-F. Li, G.-C. Guo, and A. del Campo, Commun. Phys. 3, 44 (2020).
  * Li  _et al._ (2023) B.-W. Li, Y.-K. Wu, Q.-X. Mei, R. Yao, W.-Q. Lian, M.-L. Cai, Y. Wang, B.-X. Qi, L. Yao, L. He, Z.-C. Zhou, and L.-M. Duan, PRX Quantum 4, 010302 (2023).
  * 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. Vuletić, and M. D. Lukin, Nature 568, 207 (2019).
  * 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. Vuletić, and M. D. Lukin, Nature 595, 227 (2021).
  * Gardas  _et al._ (2018) B. Gardas, J. Dziarmaga, W. H. Zurek, and M. Zwolak, Sci. Rep. 8, 4539 (2018).
  * Weinberg  _et al._ (2020) P. Weinberg, M. Tylutki, J. M. Rönkkö, J. Westerholm, J. A. Åström, P. Manninen, P. Törmä, and A. W. Sandvik, Phys. Rev. Lett. 124, 090502 (2020).
  * Bando  _et al._ (2020) Y. Bando, Y. Susa, H. Oshiyama, N. Shibata, M. Ohzeki, F. J. Gómez-Ruiz, D. A. Lidar, S. Suzuki, A. del Campo, and H. Nishimori, Phys. Rev. Research 2, 033369 (2020).
  * King  _et al._ (2022) A. D. King, S. Suzuki, J. Raymond, A. Zucca, T. Lanting, F. Altomare, A. J. Berkley, S. Ejtemaee, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, J. D. Whittaker, J. Yao, R. Harris, D. A. Lidar, H. Nishimori, and M. H. Amin, Nat. Phys. 18, 1324 (2022).
  * Cincio  _et al._ (2007) L. Cincio, J. Dziarmaga, M. M. Rams, and W. H. Zurek, Phys. Rev. A 75, 052321 (2007).
  * del Campo (2018) A. del Campo, Phys. Rev. Lett. 121, 200601 (2018).
  * del Campo  _et al._ (2022) A. del Campo, F. J. Gómez-Ruiz, and H.-Q. Zhang, Phys. Rev. B 106, L140101 (2022).
  * Gómez-Ruiz  _et al._ (2020) F. J. Gómez-Ruiz, J. J. Mayo, and A. del Campo, Phys. Rev. Lett. 124, 240602 (2020).
  * Mayo  _et al._ (2021) J. J. Mayo, Z. Fan, G.-W. Chern, and A. del Campo, Phys. Rev. Research 3, 033150 (2021).
  * del Campo  _et al._ (2021) A. del Campo, F. J. Gómez-Ruiz, Z.-H. Li, C.-Y. Xia, H.-B. Zeng, and H.-Q. Zhang, JHEP 2021, 61 (2021).
  * Gómez-Ruiz  _et al._ (2022) F. J. Gómez-Ruiz, D. Subires, and A. del Campo, Phys. Rev. B 106, 134302 (2022).
  * Ellis (2006) R. Ellis,  _Entropy, Large Deviations, and Statistical Mechanics_ , Classics in Mathematics (Springer, 2006).
  * Touchette (2009) H. Touchette, Phys. Rep. 478, 1 (2009).
  * Gambassi and Silva (2012) A. Gambassi and A. Silva, Phys. Rev. Lett. 109, 250602 (2012).
  * Goold  _et al._ (2018) J. Goold, F. Plastina, A. Gambassi, and A. Silva, “The role of quantum work statistics in many-body physics,” in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, edited by F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso (Springer International Publishing, Cham, 2018) pp. 317–336.
  * Perfetto  _et al._ (2019) G. Perfetto, L. Piroli, and A. Gambassi, Phys. Rev. E 100, 032114 (2019).
  * (45) See Supplementary Information.
  * Damski and Rams (2013) B. Damski and M. M. Rams, J. Phys. A: Math. Theor. 47, 025303 (2013).
  * Damski and Zurek (2006) B. Damski and W. H. Zurek, Phys. Rev. A 73, 063405 (2006).
  * Vershynin (2018) R. Vershynin,  _High-Dimensional Probability: An Introduction with Applications in Data Science_, Cambridge Series in Statistical and Probabilistic Mathematics (Cambridge University Press, 2018).
  * Vodola  _et al._ (2014) D. Vodola, L. Lepori, E. Ercolessi, A. V. Gorshkov, and G. Pupillo, Phys. Rev. Lett. 113, 1 (2014).
  * Yang  _et al._ (2022) J. Yang, S. Pang, A. del Campo, and A. N. Jordan, Phys. Rev. Research 4, 013133 (2022).
  * Dutta and Dutta (2017) A. Dutta and A. Dutta, Phys. Rev. B 96, 125113 (2017).
  * Pezzè  _et al._ (2017) L. Pezzè, M. Gabbrielli, L. Lepori, and A. Smerzi, Phys. Rev. Lett. 119, 250401 (2017).

Supplementary Information

##  I Fermionic base for the Transverse-Field Quantum Ising Model

To fix the notation, we briefly review how the transverse-field quantum Ising
model (TFQIM) can be mapped onto a set of non-interacting two-level systems
Dziarmaga (2010). Let us start from the chain Hamiltonian

$H[g(t)]=-J\sum_{l=1}^{N}\left[g(t)\sigma_{l}^{x}+\sigma_{l}^{z}\sigma_{l+1}^{z}\right].$
(23)

First, by applying on each site the unitary gate

$U=U_{\mathrm{Hadamard}}\sigma^{z}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&-1\\\
1&1\end{pmatrix},$ (24)

$H$ can be brought into the equivalent form

$H[g(t)]=J\sum_{l=1}^{N}\left[g(t)\sigma_{l}^{z}-\sigma_{l}^{x}\sigma_{l+1}^{x}\right].$
(25)

Second, we apply a Jordan-Wigner transformation and we introduce the fermionic
operators

$c_{l}^{\dagger}:=e^{{\mathrm{i}}\pi\Sigma_{l}}\sigma_{l}^{+},\qquad
c_{l}:=e^{-{\mathrm{i}}\pi\Sigma_{l}}\sigma_{l}^{-},\qquad\sigma_{l}^{z}=2c_{l}^{\dagger}c_{l}^{\phantom{\dagger}}-1,$
(26)

where $\Sigma_{l}$ is the string

$\Sigma_{l}:=\sum_{j=1}^{l-1}c_{j}^{\dagger}c_{j}^{\phantom{\dagger}}.$ (27)

The Hamiltonian (25) in terms of the fermionic operators reads

$H[g(t)]=-J\sum_{l}\left[c_{l}^{\dagger}c_{l+1}^{\dagger}+c_{l}^{\dagger}c_{l+1}^{\phantom{\dagger}}-c_{l}^{\phantom{\dagger}}c_{l+1}^{\dagger}-c_{l}c_{l+1}-2g(t)c_{l}^{\dagger}c_{l}^{\phantom{\dagger}}\right]-NJg(t).$
(28)

Next, we define the Fourier basis

$\tilde{c}_{k}:=\frac{1}{\sqrt{N}}\sum_{l}e^{-{\mathrm{i}}kl}c_{l},\qquad\tilde{c}_{k}^{\dagger}:=\frac{1}{\sqrt{N}}\sum_{l}e^{{\mathrm{i}}kl}c_{l}^{\dagger},$
(29)

where the symmetric definition of the Fourier transform (i.e., using the
$1/\sqrt{N}$ prefactor) is essential to have the commutation relations
$\\{\tilde{c}_{k}^{\phantom{\dagger}},\tilde{c}_{q}^{\dagger}\\}=\delta_{kq}$
without spurious factors of $N$. One needs to impose the boundary conditions
to fix the values that $k$ can acquire Damski and Rams (2013). We assume
periodic boundary conditions in the spin chain, i.e.,
$\vec{\sigma}_{l+N}=\vec{\sigma}_{l}$. In the fermionic representation, there
is a dependence on the number of particles. If it is even, a particle must
acquire a negative phase when performing a full circle. This means
$P^{+}c_{l+N}P^{+}=-c_{l}$, where $P^{+}$ projects onto the even sector.
Similarly, $P^{-}c_{l+N}P^{-}=c_{l}$. We choose to work with an even number of
particles, whence

$c_{l+N}=-c_{l}\quad\implies\quad\frac{1}{\sqrt{N}}\sum_{k}e^{{\mathrm{i}}kl}\tilde{c}_{k}=\frac{1}{\sqrt{N}}\sum_{k}e^{{\mathrm{i}}k(l+N)}\tilde{c}_{k}\quad\implies\quad
e^{{\mathrm{i}}kN}=-1.$ (30)

One can verify that the solution is

$k=\frac{2\pi}{N}\left(n+\frac{1}{2}\right)\qquad\text{for
}n=-\frac{N}{2},-\frac{N}{2}+1,\dots,\frac{N}{2}-1.$ (31)

At this point, the fermionic Hamiltonian rewritten in momentum space reads

$H[g(t)]=2J\sum_{k>0}\left[-{\rm i}\sin
k\left(\tilde{c}_{k}^{\dagger}\tilde{c}_{-k}^{\dagger}+\tilde{c}_{-k}\tilde{c}_{k}\right)+\left(g(t)-\cos
k\right)\left(\tilde{c}_{k}^{\dagger}\tilde{c}_{k}^{\phantom{\dagger}}-\tilde{c}_{-k}^{\phantom{\dagger}}\tilde{c}_{-k}^{\dagger}\right)\right].$
(32)

Alternatively,

$H[g(t)]=2J\sum_{k>0}\begin{pmatrix}\tilde{c}_{k}^{\dagger},&\tilde{c}_{-k}^{\phantom{\dagger}}\end{pmatrix}\begin{pmatrix}g(t)-\cos
k&-{\rm i}\sin k\\\ {\rm i}\sin k&-g(t)+\cos
k\end{pmatrix}\begin{pmatrix}\tilde{c}_{k}\\\
\tilde{c}_{-k}^{\dagger}\end{pmatrix}$ (33)

and, defining

$\psi_{k}^{\dagger}:=\begin{pmatrix}\tilde{c}_{k}^{\dagger},&\tilde{c}_{-k}^{\phantom{\dagger}}\end{pmatrix},\qquad\psi_{k}:=\begin{pmatrix}\tilde{c}_{k}\\\
\tilde{c}_{-k}^{\dagger}\end{pmatrix},$ (34)

the Hamiltonian acquires the compact form

$H[g(t)]=2J\sum_{k>0}h_{k}[g(t)],$ (35)

where

$h_{k}[g(t)]\equiv\psi_{k}^{\dagger}\left[(g(t)-\cos k)\tau^{z}+\sin
k\,\tau^{y}\right]\psi_{k},$ (36)

and the Pauli matrices $\tau^{x,y,z}$ act on the two entries of the
$\psi_{k}$’s. Notice that $\psi_{k}$ has only positive $k$’s in the subscript.
Equation (35) is exactly the same as Eq. (2) in the main text.

##  II Diagonalization of the Transverse-Field Quantum Ising Model

We want to diagonalize the Hamiltonian (2) through a Bogoliubov rotation.
Thus, we change the basis for the creation/annihilation operators as

$\begin{pmatrix}\gamma_{k}(t)\\\
\gamma_{-k}^{\dagger}(t)\end{pmatrix}=\begin{pmatrix}\cos[\theta_{k}(t)/2]&-{\mathrm{i}}\sin[\theta_{k}(t)/2]\\\
-{\mathrm{i}}\sin[\theta_{k}(t)/2]&\cos[\theta_{k}(t)/2]\end{pmatrix}\begin{pmatrix}\tilde{c}_{k}\\\
\tilde{c}_{-k}^{\dagger}\end{pmatrix}.$ (37)

Again, only $k>0$ are used, leaving an explicit minus sign where needed. The
value of $\theta_{k}$ is fixed by the diagonalization of $h_{k}$, which yields
the relation

$\cos\theta_{k}(t)=\frac{g(t)-\cos k}{\sqrt{g^{2}(t)+1-2g(t)\cos k}}.$ (38)

It follows that

$H[g(t)]=\sum_{k>0}\begin{pmatrix}\gamma_{k}^{\dagger}(t),&\gamma_{-k}(t)\end{pmatrix}\begin{pmatrix}\epsilon_{k}(t)&0\\\
0&-\epsilon_{k}\end{pmatrix}\begin{pmatrix}\gamma_{k}(t)\\\
\gamma_{-k}^{\dagger}(t)\end{pmatrix},$ (39)

with

$\epsilon_{k}(t)=2J\sqrt{g^{2}(t)+1-2g(t)\cos k}.$ (40)

Using $\\{\gamma,\gamma^{\dagger}\\}=1$, one can rewrite the Hamiltonian as

$H[g(t)]=\sum_{k}\epsilon_{k}(t)\left(\gamma^{\dagger}_{k}(t)\gamma^{\phantom{\dagger}}_{k}(t)-\frac{1}{2}\right),$
(41)

where the sum runs over both positive and negative $k$’s.

##  III Additional plots

We perform additional numerical analysis of the kink statistics in the driven
TFQIM, establishing its regime of universality. In Fig. 2(a) we plot the
average density of kink pairs. As predicted by the conventional KZM, the
power-law scaling $\langle\hat{\rho}_{N}\rangle\propto\tau_{Q}^{-1/2}$ holds
if quench times are not too long, i.e., before the onset of adiabaticity. From
this figure, it can be checked that the parameter values used to make Fig. 1
lie well within the KZ scaling region, characterized by the universal power-
law scaling with the quench time.

Figure 2: (a) Average kink-pair density produced by the finite-time crossing
of the phase transition. One can see the curves collapsing on the single
power-law scaling $\langle\hat{\rho}_{N}\rangle\sim\tau_{Q}^{-1/2}$ as the
quench time is increased. However, at finite $N$ and for quench times too
large, the KZM breaks down since the dynamics becomes purely adiabatic. (b)
Cumulants of the probability distribution $P(\rho)$. As already observed in
Refs. del Campo (2018); Mayo  _et al._ (2021); Cui  _et al._ (2020); Bando
_et al._ (2020); King  _et al._ (2022), the cumulants scale with the same
power law as a function of the quench time, a hallmark of the underlying
Poisson-binomial distribution. (c) The ratios between cumulants fastly
approach the universal values $\kappa_{2}/\kappa_{1}=0.293\dots$,
$\kappa_{3}/\kappa_{1}=0.0334\dots$, $\kappa_{3}/\kappa_{2}=0.114\dots$, shown
as dotted lines. To make this figure, the same quench protocol of Fig. 1 was
employed.

To further establish the regime governed by a universal power-law scaling, in
Fig. 2(b)–(c) we compare the scaling of the first three cumulants of the kink-
pair probability distribution $P(\rho)$, as a function of the quench time
$\tau_{Q}$. Making use of the fact that the kink-pair distribution is Poisson-
Binomial, the explicit expression for the first few cumulants in the TFQIM
reads del Campo (2018); Cui  _et al._ (2020); Bando  _et al._ (2020)

$\kappa_{1}=\sum_{k>0}p_{k},\qquad\kappa_{2}=\sum_{k>0}p_{k}(1-p_{k}),\qquad\kappa_{3}=\sum_{k>0}p_{k}(1-p_{k})(1-2p_{k}).$
(42)

We obtain the excitation probabilities $p_{k}$ numerically without relying on
the LZ approximation. The universal full counting statistics of topological
defects beyond KZM del Campo (2018); Mayo  _et al._ (2021); Cui  _et al._
(2020); Bando  _et al._ (2020); King  _et al._ (2022) dictates that all
cumulants exhibit a universal power-law scaling with the quench rate, with
cumulant ratios being constant. Deviations from this universal behavior are
apparent for moderate quench rate satisfying $J\tau_{Q}/\hbar<10$, with high-
order cumulants being more sensitive than the mean to nonuniversal effects.
Panels (b)–(c) of Fig. 2 identify the regime in which the universal power-law
scaling holds, not only for the mean but also for higher-order cumulants.

Figure 3: The rate function $\bar{I}$ can be equivalently accessed by
computing directly the logarithm of the probability density function. However,
even for moderately large $N$, the discreteness of the number of kink pairs
results in few data points. To make this figure, the same quench protocol of
Fig. 1 was employed.

We complete this supplementary analysis with Fig. 3, where we explore the
convergence of the rate function in LDT from numerically-exact simulations
using a finite system size. Specifically, we compared the numerically-
evaluated expression $-\frac{1}{N}\ln P(\bar{\rho})$ with the theoretical
asymptotic value derived in the main text. Given the finite system size $N$,
the random variable $\hat{\rho}_{N}$ can assume only the discrete values
$1/N,2/N,\dots,1/2$. Thus, only a few data points fall into the relevant
region of $\bar{\rho}$, when the data is appropriately scaled. By contrast,
the computation of $\bar{I}(\bar{\rho})$ makes use of continuous values of
$\bar{\rho}$, leading to the better collapse of Fig. 1. As a result, the rate
function obtained by numerically-exact simulations converges to the
thermodynamic limit for moderate system sizes $N=100-1000$. Deviations from
the analytical expression derived in the main text simply arise from the fact
that the probabilities $p_{k}$ deviate from the LZ approximation for moderate
quench times. This is consistent with the fact that the LZ and the validity of
the KZM power-law scaling are restricted to the limit of slow quenches.

##  IV Universality of LDT beyond KZM in systems with fast-decaying long-range
interactions

In this section, we first argue that near the critical region $g=1$, the
spectrum of excitations in the long-range deformation of the TFQIM model,
i.e., Eq. (43) is the same as the one in the TFQIM for $\alpha\geq 2$. Then we
show that the effect of the long-range is to renormalize the quench time by a
positive constant. Therefore in the slow quench limit, the universal KZM
scaling and the rate function in the LDT are robust against fast-decaying
long-range interactions. As such, we expect the rate function Eq. (21) in the
main text also holds for systems with fast decaying long-range interactions.

###  IV.1 The long-range Ising and Kitaev model

To explore the robustness of our predictions in long-range systems, we
consider a long-range Hamiltonian

$\displaystyle H_{\text{LR}}[g(t)]$ $\displaystyle=$
$\displaystyle-J\sum_{j=1}^{N}\sigma_{j}^{z}\sigma_{j+1}^{z}-Jg(t)\sum_{j=1}^{N}\sigma_{j}^{x}+\frac{J}{2}\sum_{j=1}^{N}\sum_{l=2}^{N-1}\kappa_{l,\alpha}(\sigma_{j}^{z}\sigma_{j+l}^{z}-\sigma_{j}^{y}\sigma_{j+l}^{y})\bigotimes_{n=1}^{l-1}\sigma_{j+n}^{x},$
(43)

where

$\kappa_{l,{\alpha}}=\begin{cases}\frac{1}{l^{\alpha}}&1\leq l\leq N/2\\\
\frac{1}{(N-l)^{\alpha}}&\text{otherwise},\end{cases}$ (44)

As for the short-range case, we use periodic boundary conditions
$\sigma_{1}^{a}=\sigma_{N+1}^{a}$. We note that $H_{\text{LR}}$, like its
short-ranged counterpart, possesses a $\mathbb{Z}_{2}$ symmetry, i.e.,

$\left(\bigotimes_{i}\sigma_{i}^{x}\right)H_{\text{LR}}\left(\bigotimes_{i}\sigma_{i}^{x}\right)=H_{\text{LR}}.$
(45)

As before, upon using the Hadamard gate to rotate the axes

$\sigma_{j}^{x}\to-\sigma_{j}^{z},\qquad\sigma_{j}^{y}\to\sigma_{j}^{y},\qquad\sigma_{j}^{z}\to\sigma_{j}^{x},$
(46)

one obtains

$H_{\text{LR}}[g(t)]=-J\sum_{j=1}^{N}\sigma_{j}^{x}\sigma_{j+1}^{x}+Jg(t)\sum_{j=1}^{N}\sigma_{j}^{z}+\frac{J}{2}\sum_{j=1}^{N}\sum_{l=2}^{N-1}\kappa_{l,\alpha}(\sigma_{j}^{x}\sigma_{j+l}^{x}-\sigma_{j}^{y}\sigma_{j+l}^{y})\bigotimes_{n=1}^{l-1}(-\sigma_{j+n}^{z}).$
(47)

Restricting to the even fermionic parity sector, we perform the Jordan-Wigner
transformation defined in Eq. (26) and obtain

$H_{\text{LRK}}[g(t)]=-J\sum_{j=1}^{N}(c_{j}^{\dagger}c_{j+1}+c_{j+1}^{\dagger}c_{j})+2Jg(t)\sum_{j=1}^{N}\bigg{(}c_{j}^{\dagger}c_{j}-\frac{1}{2}\bigg{)}-J\sum_{j=1}^{N-1}\sum_{l=1}^{N-j}\kappa_{l,\alpha}(c_{j}c_{j+l}+c_{j+l}^{\dagger}c_{j}^{\dagger}).$
(48)

As explained above, the boundary conditions at the level of the fermionic
Hamiltonian become anti-periodic: $c_{j}=-c_{N+j}$. Eq. (48) is the one-
dimensional long-range Kitaev (LRK) chain with power-law decaying
superconducting pairing (Vodola  _et al._ , 2014). After a Fourier transform,
Eq. (48) can be rewritten as

$H[g(t)]=2J\sum_{k>0}h_{k}[g(t)],$ (49)

where

$\displaystyle h_{k}[g(t)]$
$\displaystyle\equiv\psi_{k}^{\dagger}\left[(g(t)-\cos k)\tau^{z}+f_{\alpha
k}\tau^{y}\right]\psi_{k},$ $\displaystyle=[g(t)-\cos
k](\tilde{c}_{k}^{\dagger}\tilde{c}_{k}-\tilde{c}_{-k}\tilde{c}_{-k}^{\dagger})-\mathrm{i}f_{\alpha
k}(\tilde{c}_{k}\tilde{c}_{-k}-\tilde{c}_{-k}^{\dagger}\tilde{c}_{k}^{\dagger}),$
(50)

and the function $f_{\alpha k}$ is defined as

$f_{\alpha k}\equiv\frac{1}{2}\sum_{l=1}^{N-1}\kappa_{l,\alpha}\sin(kl),$ (51)

with the quasi-momentum taking values
$k=\pm\frac{1}{2}\frac{2\pi}{N},\,\pm\frac{3}{2}\frac{2\pi}{N}\cdots,\,\pm\frac{2\pi}{N}(\frac{N}{2}-\frac{1}{2})$.
In the case of $\alpha=\infty$, it holds that $f_{\infty k}=\sin k$, and the
model reduces to the short-range Kitaev chain, i.e., the TFQIM after Jordan-
Wigner-transforming back.

Equation (49) can be further diagonalized through the Bogoliubov
transformation in the form

$H[g(t)]=\sum_{k>0}\epsilon_{k}(t)\psi_{k}^{\dagger}U_{k}^{\dagger}(t)\tau^{z}U_{k}(t)\psi_{k},$
(52)

where

$\epsilon_{k}(t)=2J\sqrt{[f_{\alpha k}]^{2}+[g(t)-\cos k]^{2}},$ (53)

and

$U_{k}(t)\equiv\begin{pmatrix}\cos\left[\theta_{k}(t)/2\right]&-\text{i}\sin\left[\theta_{k}(t)/2\right]\\\
-\text{i}\sin\left[\theta_{k}(t)/2\right]&\cos\left[\theta_{k}(t)/2\right]\end{pmatrix},$
(54)

with

$\sin\theta_{k}(t)=\frac{f_{\alpha k}}{\sqrt{[f_{\alpha k}]^{2}+[g(t)-\cos
k]^{2}}},\quad\cos\theta_{k}(t)=\frac{[g(t)-\cos k]}{\sqrt{[f_{\alpha
k}]^{2}+[g(t)-\cos k]^{2}}}.$ (55)

Denoting

$\begin{pmatrix}\gamma_{k}(t)\\\ \gamma_{-k}^{\dagger}(t)\end{pmatrix}\equiv
U_{k}(t)\begin{pmatrix}\tilde{c}_{k}\\\
\tilde{c}_{-k}^{\dagger}\end{pmatrix},$ (56)

the Hamiltonian can be rewritten as

$H[g(t)]=\sum_{k}\epsilon_{k}(t)\left[\gamma_{k}^{\dagger}(t)\gamma_{k}(t)-\frac{1}{2}\right],$
(57)

where the summation is over all the modes, both positive and negative.

###  IV.2 Properties of $f_{\alpha k}$

Let us list a few properties of $f_{\alpha k}$. First of all, $f_{\alpha k}$
is odd in $k$, i.e.,

$f_{\alpha k}=-f_{\alpha,-k}.$ (58)

Note that in the thermodynamic limit, the quasi-momentum $k\in[-\pi,\,\pi]$.
While $\lim_{k\to\pm\pi}f_{\alpha k}=0$, one should note that $\lim_{k\to
0}f_{\alpha}(k)$ does not necessarily vanish. In general, if
$\kappa_{l,\alpha}$ decays slowly as $l$ increases, $f_{\alpha k}$ may become
singular at $k=0$, as shown in Ref. (Yang  _et al._ , 2022). In fact, one can
express $f_{\alpha k}$ in terms of the polylogarithm function

$f_{\alpha k}=-{\rm i}\left[\text{Li}_{\alpha}(e^{{\rm
i}k})-\text{Li}_{\alpha}(e^{-{\rm i}k})\right].$ (59)

Using the expansion of polylogarithm function, one can show that
$\alpha\notin\mathbb{N}^{+}$, $f_{\alpha k}\propto 1/k^{1-\alpha}$ (Vodola
_et al._ , 2014), which is adopted by Ref. (Dutta and Dutta, 2017) to discuss
KZM in the long-range Kitaev chain. We note that this scaling behavior can be
alternatively proved using the Euler-Maclaurin formula for any $\alpha\geq 0$
without the constraint $\alpha\notin\mathbb{N}^{+}$. Thus, in view of Eq.
(58), $f_{\alpha k}$ allows the following expansion near $k=0$,

$f_{\alpha k}=\frac{\xi_{0}}{k^{1-\alpha}}+\xi_{1}k+O(k^{3}).$ (60)

###  IV.3 LDT for the kink density with fast-decaying long-range interactions

The phase diagram of the long-range Kitaev chain is very rich (Vodola  _et
al._ , 2014; Pezzè  _et al._ , 2017), including in particular topological
phase transitions. With the spectrum Eq. (53), for $\alpha>1$ one can easily
verify that the gap closes only when $g=\pm 1$, with the gapless mode being
$k=0$. Furthermore, according to Eq. (60), if one only considers $\alpha\geq
2$, then

$f_{\alpha k}=k+\mathcal{O}(k^{3}).$ (61)

This is the same behaviour of the spectrum of the TFQIM near $k=0$, where
$f_{\alpha k}=\sin k$. Clearly, for both $-1<g<1$ and $g>1$, one can
continuously deform $\alpha$ from $2$ to $\infty$ without closing the gap. As
we will show below, the kink density statistics and the large deviations
theory are the same in the slow quench limit. We consider the ramp in the main
text, i.e.,

$g(t)=g_{c}\left(1-\frac{t}{\tau_{Q}}\right),$ (62)

starting at an initial time before $-\tau_{Q}$, outside the freeze-out region
del Campo and Zurek , and ending at $\tau_{Q}$, whence $g(\tau_{Q})=0$. We
investigate the statistics of kink density and work with the adiabatic basis
at $t=-\infty$,

$\displaystyle|0\rangle,\qquad|k,-k\rangle=c_{k}^{\dagger}c_{-k}^{\dagger}|0\rangle=\gamma_{k}^{\dagger}(-\infty)\gamma_{-k}(-\infty)|0\rangle,$
(63)
$\displaystyle|k,0\rangle=c_{k}^{\dagger}|0\rangle=\gamma_{k}^{\dagger}(-\infty)|0\rangle,\qquad|0,-k\rangle=c_{-k}^{\dagger}|0\rangle=\gamma_{-k}(-\infty)|0\rangle.$
(64)

For each mode Hamiltonian $Jh_{k}[g(t)]$ preserves the fermionic parity, and
excitations remain in the subspace spanned by $|0\rangle,|k,-k\rangle$.
Projecting Eq. (50) onto this subspace, one obtains the effective Hamiltonian

$Jh_{k}[g(t)]=J\begin{pmatrix}[g(t)-\cos k]&-{\mathrm{i}}f_{\alpha k}\\\
{\mathrm{i}}f_{\alpha k}&-[g(t)-\cos k]\end{pmatrix}.$ (65)

The excitation generated during the process can be characterized by the
Landau-Zener formula, provided $|k|<\pi/2$:

$p_{k}=\left\langle\gamma_{k}^{\dagger}(\tau_{Q})\gamma_{k}(\tau_{Q})\right\rangle\approx\exp\bigg{(}-\frac{2\pi
J\tau_{Q}}{\hbar}f_{\alpha}^{2}(k)\bigg{)},$ (66)

where the average is over the state of the system at the final time
$t=\tau_{Q}$. The scaled cumulant generating function becomes resembles that
in the TFQIM

$\lambda(\theta)=\frac{1}{2\pi}\int_{0}^{\pi}dk\ln\bigg{[}1+(e^{\theta}-1)\exp\bigg{(}-\frac{2\pi
J\tau_{Q}}{\hbar}f_{\alpha}^{2}(k)\bigg{)}\bigg{]},$ (67)

where we have extended the integral up to $\pi$, an approximation with
exponential accuracy. When $\alpha\geq 2$, according to Eq. (60),

$f_{\alpha k}=\xi k+O(k^{3}).$ (68)

In the slow quench limit where $\tau_{Q}\to\infty$, one can perform a Gaussian
approximation, leading to

$p_{k}=\left\langle\gamma_{k}^{\dagger}(\tau_{Q})\gamma_{k}(\tau_{Q})\right\rangle=\exp\bigg{(}\frac{-2\pi
J\tau_{Q}\xi^{2}k^{2}}{\hbar}\bigg{)},$ (69)

and thus

$\lambda(\theta)=\frac{1}{2\pi}\int_{0}^{\pi}dk\ln\bigg{[}1+(e^{\theta}-1)\exp\bigg{(}\frac{-2\pi
J\tau_{Q}\xi^{2}k^{2}}{\hbar}\bigg{)}\bigg{]}.$ (70)

Comparing the results for TFQIM, we see that the effect of the fast decaying
long-range interactions is to renormalize the quench time by a positive
constant, keeping the structure of the scaled cumulant generating function
otherwise unchanged.