Kibble-Zurek mechanism in the Ising Field Theory

The Kibble-Zurek mechanism captures universality when a system is driven through a continuous phase transition. Here we study the dynamical aspect of quantum phase transitions in the Ising Field Theory where the critical point can be crossed in different directions in the two-dimensional coupling space leading to different scaling laws. Using the Truncated Conformal Space Approach, we investigate the microscopic details of the Kibble-Zurek mechanism in a genuinely interacting field theory. We demonstrate dynamical scaling in the non-adiabatic time window and provide analytic and numerical evidence for specific scaling properties of various quantities, including higher cumulants of the excess heat.


Introduction
The Kibble-Zurek mechanism (KZM) describes the dynamical aspects of phase transitions and captures the universal features of nonequilibrium dynamics when a system is driven slowly across a continuous phase transition. The original idea is due to Kibble, who studied cosmological phase transitions in the early Universe [1,2]. He showed that as the Universe cools below the symmetry breaking temperature, instead of perfect ordering, domains form whose typical size depends on the cooling rate. Later Zurek pointed out that this phenomenon can be observed in condensed matter systems as well [3,4]. The physical mechanism originates in the fact that at a critical point both the correlation length and the correlation time (relaxation time) diverge, leading to an inevitable breakdown of adiabaticity. As a consequence, the final state will not be perfectly ordered but will consist of domains with different symmetry breaking orders separated by defects or domain walls. However, in the process a typical time scale and a corresponding length scale emerges related to the instant when the system deviates from the adiabatic course. These quantities, diverging as the rate at which the phase transition is crossed approaches zero, are the only scales in the problem. As a consequence, the density of domain walls as well as other quantities obey scaling laws in terms of the speed of the ramp.
It is a natural question whether the same phenomena occur also at zero temperature, i.e. for quantum critical points. A systematic study of the KZM in quantum phase transitions started with the works [5,6,7,8]. Quantum phase transitions are different from transitions at finite temperature: they correspond to a qualitative change in the ground state of a quantum system and are driven by quantum fluctuations. Importantly, the time evolution is unitary and there is no dissipation. In spite of these differences, the scaling behaviour essentially coincides with the classical case [5,6,7,8,9,10]. The scaling behaviour was extended to other observables beyond the defect density to correlation functions [11,12,13], entanglement entropy [13,14,15], excess heat [16,17,18], and also to different ramp protocols [10,16,19], including quenches from the ordered to the disordered phase. The scaling laws can also be derived using the framework of adiabatic perturbation theory [7,20,21,16,17,22,19,23]. The reader interested in the KZM in the context of quantum phase transitions is referred to the excellent reviews [24,25,26].
The simplest approximation which leads to the right scaling exponents assumes that when adiabaticity is lost, the system becomes completely frozen and reenters the dynamics only some time after crossing the critical point. This freeze-out scenario or impulse approximation has been refined recently by taking into account the actual evolution of the system in the non-adiabatic time window [27,28,29,30,31,15,32,33]. Since the Kibble-Zurek length and time scales are the only relevant scales, the non-adiabatic evolution features dynamical scaling, i.e. the time dependence of various observables is given by scaling functions.
The Kibble-Zurek mechanism was also extended beyond the mean values to the full statistics of observables. The number distribution of defects was computed in the Ising chain [13,34] and was argued to exhibit universality [35]. Similarly, the work statistics and its cumulants were also studied and found to satisfy scaling relations [36,37,38].
In this work we aim to study different aspects the quantum Kibble-Zurek mechanism in a simple but nontrivial field theory, the paradigmatic Ising Field Theory. This theory is an ideal testing ground as it allows one to study both free and genuinely interacting integrable systems. Our motivation for this choice is twofold. First, we wish to study the KZM in a field theory at the microscopic level of states. Second, we would like to test the recent predictions for the universal dynamical scaling and the scaling behaviour of the higher cumulants of the work in an interacting model.
As we focus on an interacting theory, we need to use a numerical tool for our studies. We use a nonperturbative numerical method, the Truncated Space Approach [63,64,65]. Apart from its long-standing history to capture equilibrium properties of perturbed conformal field theories [66,67,68,69,70,71,72,73,74,75,76,77,78], recent applications demonstrate that it is capable to describe non-equilibrium dynamics in different models [79,80,81,82,83,84]. This approach gives us access to the microscopic data and full statistics of observables so we can investigate the KZM at work at the lowest level, and being nonperturbative and independent of integrability, it allows us to study the dynamics of the interacting field theory.
The paper is organised as follows. In Sec. 2 we outline the context of our work and review the scaling laws predicted by the Kibble-Zurek mechanism for quantum phase transitions. We proceed by defining the model in which we study the Kibble-Zurek mechanism and discuss the adiabatic perturbation theory that provides another viewpoint on the scaling laws. The main body of the text presents an in-depth analysis of the Kibble-Zurek mechanism in the Ising Field Theory. In Sec. 3 we explore the implications of driving a system across a critical point on the statistics of work function and examine the behaviour of energy eigenstates to check the hypothesis of the KZM at a fundamental level. Sec. 4 discusses the dynamical critical scaling with the time and length scales corresponding to the deviation from the adiabatic course and demonstrates that the KZ scaling can be observed in the interacting E 8 model. In Sec. 5 we show that the appearance of the scaling connected to the Kibble-Zurek mechanism is not limited to local observables but it is present also in higher cumulants of the distribution of the excess heat. Finally, Sec. 6 finishes the paper with concluding remarks and possible future directions. Technical details concerning the relation of the adiabatic perturbation theory to the E 8 model, the scaling limit of the analytic solution of the dynamics on the transverse field Ising chain and the applicability of TCSA to the study of KZM are discussed in the Appendices.

Model and methods
In this section we describe the context of our work by introducing the concepts of the universal non-adiabatic behaviour that manifests itself in power-law dependence of several quantities on the time scale of the non-equilibrium ramp protocol, known under the name of Kibble-Zurek scaling. Then we discuss the model in which we study the KZ scaling, the Ising Field Theory which is the low energy effective theory of the transverse field Ising chain in the vicinity of its critical point. After introducing its main properties, we address the methods that are going to be used to examine the Kibble-Zurek scaling. In the limit of slow ramps, one can employ a perturbative approach, the adiabatic perturbation theory (APT) to investigate the time evolution. We give an overview of this approach, focusing on its application to universal dynamics near quantum critical points. The non-equilibrium dynamics of the Ising Field Theory is amenable to an efficient numerical non-perturbative treatment based on the truncated conformal space approach (TCSA), which we review briefly at the end of the section.

The Kibble-Zurek mechanism
In this section we summarise the KZ scaling laws in a fairly general fashion. Let us consider a perturbation of a quantum critical point (QCP) by some operator with scaling dimension ∆. The strength of the perturbation is characterised by a coupling constant δ with δ = 0 corresponding to the critical point. Imagine that we prepare the system in its ground state and drive it through its QCP by changing δ in time, i.e. by performing a ramp. For the sake of generality, we consider ramps that cross the phase transition in a power-like fashion, i.e. near the QCP where τ Q is the rate of the quench. The essence of the KZM is that due to the divergence of the relaxation time of the system at the QCP, known as critical slowing down, the system cannot follow adiabatically the change no matter how slow it is, and falls out of equilibrium meaning that it will be in an excited state with respect to the instantaneous Hamiltonian. However, due to universality near the critical point the time and length scales corresponding to the deviation from the adiabatic course depend on the quench rate τ Q as a power-law. The scaling can be determined by the following simple argument. The correlation length diverges in the phase transition corresponding to this particular perturbation as ξ ∝ δ −ν where ν is the standard equilibrium critical exponent related to the scaling dimension ∆ of the perturbing operator by ν = (2 − ∆) −1 . Similarly, the correlation or relaxation time diverges as ξ t ∝ ξ z ∝ δ −νz , where z is the dynamical critical exponent. If the change of ξ t within a relaxation time is much smaller than the relaxation time itself,ξ t ξ t ξ t , then the evolution is adiabatic. This is the case for times However, once we reach t ≈ −τ KZ , the rate of change of the correlation time becomesξ t ≈ 1 and the evolution becomes non-adiabatic. At this Kibble-Zurek time τ KZ , the correlation time scales with the quench rate τ Q as τ KZ itself: 3) The first formulation of Kibble-Zurek arguments depicted the non-adiabatic interval of time evolution as a simple freeze-out referring to the assumption that the state is literally frozen in the non-adiabatic regime t ∈ [−τ KZ , τ KZ ]. At t = τ KZ on the other side of the QCP, the system finds itself in an excited state with correlation length ξ KZ = ξ(−τ KZ ). If the system is now in the ordered phase, it implies that the typical linear size of the ordered domains are ∼ ξ KZ , so the density of excitations corresponding to defects (domain walls) in spatial dimension d is Recently, the freeze-out scenario was refined by taking into account the evolution of the system and change of the correlation length in the time interval −τ KZ < t < τ KZ [27,28,29,31]. The latter is caused by moving domain walls at the typical velocity corresponding to their typical wave number k ∼ ξ −1 KZ and energy ε(k) ∼ k z ∼ ξ −z KZ . The velocity of this "sonic horizon" [31] is The correlation length at t = τ KZ is then which, due to Eq. (2.3), is proportional to ξ KZ . This means that ξ KZ is still the only relevant length scale so the scaling laws are not altered. Still, nontrivial predictions can be made concerning the non-adiabatic or impulse region −τ KZ < t < τ KZ [29,31,32] due to the fact that the KZ time and correlation length, τ KZ and ξ KZ , are the only relevant scales for a slow enough ramp protocol. Consequently, timedependent correlation functions are described in terms of scaling functions of the rescaled variables t/τ KZ and x/ξ KZ in the KZ scaling limit τ KZ → ∞. For example, one-and twopoint functions of an operator (2.6) where F and G are scaling functions depending on the operator O and we assumed translational invariance. Note that for one-point functions the scaling holds in the adiabatic regime t < −τ KZ as well, since there the expectation value depends only on the distance from the critical point, which is the function of the dimensionless time t/τ Q : where in the last step we used the relation (2.2).
Considering the generic nature of arguments presented above it is tempting to ask how precisely they describe the actual non-equilibrium dynamics of quantum systems. The scaling relations are supported by exact calculations in the free fermionic Ising chain where the dynamics of low-energy modes can be mapped to the famous Landau-Zener transition problem [85,5,8,31]. In other quantum phase transitions, when exact solutions are not available, the scaling can be analysed by a perturbative expansion in the derivative of the time-dependent coupling as a small parameter. This approach that uses adiabatic perturbation theory predicts the same scaling as the arguments of Kibble-Zurek mechanism in several models besides the Ising chain [7,21,17,19]. This formalism is useful to apply the generic scaling arguments outside the non-adiabatic regime for quantities that are beyond the scope of the initial formulation of KZM [38]. Together with the non-perturbative numerical method employed in our work it can be used to establish the validity of the scaling relations listed above for an interacting model as well.
To do so, we have to address the question of finite size effects. These are of importance due to the fact that the TCSA method requires finite volume, while the arguments presented above make use of a divergent length scale ξ KZ . Clearly, finite volume can bring about adiabatic behaviour if where ξ and ξ t are the correlation length and time at the initial state. If the quench rate τ Q is significantly larger than this, the transition is adiabatic due to the fact that finite volume opens the gap. One way to compensate this effect is the rescaling of the volume parameter with the appropriate power of the quench rate [28]. However, if then the finite size effects are negligible. As we are going to illustrate in Sec. 3.3, we can set the parameters of the numerical TCSA method such that this relation is satisfied and there is no need to rescale the volume parameter.

KZM in the Ising Field Theory
After setting up the context of our work, we now turn to the model in consideration: the Ising Field Theory that is the scaling limit of the critical transverse field Ising chain. The Hamiltonian of the latter reads where σ α i with α = x, y, z are the Pauli matrices at site i, the strength of the ferromagnetic coupling J sets the energy scale, and h x J and h z J are the longitudinal and transverse magnetic fields, respectively. We set periodic boundary conditions, σ α L+1 = σ α 1 . The model is fully solvable in the absence of the longitudinal field, h x = 0, when it can be mapped to free Majorana fermions via the nonlocal Jordan-Wigner transformation. The Hilbert space is composed of two sectors based on the conserved parity of the fermion number. The fermionic Hamiltonian will be local provided we impose anti-periodic boundary conditions for the fermionic operators in the even Neveu-Schwarz (NS) sector and periodic boundary conditions in the odd Ramond (R) sector.
The transverse field Ising model is a paradigm of quantum phase transitions: in infinite volume, for h z < 1 the ground state manifold is doubly degenerate, spontaneous symmetry breaking selects the states (|0 NS ± |0 R )/ √ 2 with finite magnetisation σ = ±(1 − h 2 z ) 1/8 (here |0 NS/R are the ground states in the two sectors). In finite volume, there is an energy split between the states |0 NS and |0 R which is exponentially small in the volume, and the ground state is |0 NS . In the paramagnetic phase for h z > 1, the ground state is always |0 NS and the magnetisation vanishes. The quantum critical point (QCP) separating the ordered and disordered phases is located at h z = 1, which can also be seen from the behaviour of the gap, ∆ = 2J|1 − h z |, vanishing at the QCP. In the ferromagnetic phase, the massive fermionic excitations can be thought of domain walls separating domains of opposite magnetisations, and with periodic boundary conditions their number is always even 1 . In the paramagnetic phase the excitations are essentially spin flips in the z direction.
For h x = 0 the model is not integrable 2 for any value of h z , but features weak confinement: the nonzero longitudinal field splits the degeneracy between the two ground states with an energy difference proportional to the system size. The domain walls cease to be freely propagating excitations, as the energy cost increases with the distance between two neighbouring domain walls that have a domain of the wrong magnetisation between them. The new excitations are a tower of bound states, sometimes called 'mesons' in analogy with quark confinement in the strong interaction.
The low energy effective theory describing the model near the critical point is the Ising field theory, obtained in the scaling limit J → ∞, a → 0, h z → 1 such that speed of light c = 2Ja and the gap ∆ = 2J|1 − h z | are fixed (a is the lattice spacing). The critical point is described by a conformal field theory (CFT) of free massless Majorana fermions having central charge c = 1/2. Due to relativistic invariance, the dynamical critical exponent is z = 1. The two relevant operators in this CFT are the magnetisation σ (scaling dimension 1/8) and the 'energy density' ε (scaling dimension 1), giving rise to the two relevant perturbations corresponding to the magnetic fields of the lattice model. The Hamiltonian of the resulting field theory finite volume L is given by The precise relations between the lattice and continuum versions of the magnetic field and the magnetisation operator are For h = 0 the Hamiltonian describes the dynamics of a free Majorana fermionic field with mass |M | (we set the speed of light to one, c = 1). We will refer to this choice of parameters in the M − h parameter plane of the theory (2.11) as the "free fermion line" (see Fig. 2.1). The QCP at M = 0 separates the paramagnetic phase M > 0 from the ferromagnetic phase M < 0. The coupling is proportional to the mass gap and since the correlation length is the inverse of the gap, ν = 1.
Interestingly, there is another set of parameters that corresponds to an integrable field theory: M = 0 with h finite 3 . The spectrum of this theory can be described in terms of eight stable particles, the mass ratios and scattering matrices of which can be written in terms of the representations of the exceptional E 8 Lie group. From now on, we are going to refer to this specific set of parameters as the "E 8 integrable line" (see Fig. 2.1). The lightest particle with mass m 1 sets the energy scale which is connected to the coupling h as m 1 = (4.40490857 . . . )|h| 8/15 . (2.14) The exponent reflects that along the E 8 line (σ perturbation) ν = 8/15 and z = 1. Moving particle states are built up as combinations of particles with finite momenta from the same or different species.
In the following we are going to consider ramp protocols along the integrable lines, indicated by the coloured arrows in Fig. 2.1, where one of the couplings is varied such that the system crosses the critical point at a constant rate, corresponding to a linear ramp profile, where λ stands for M or h and the other coupling is set to zero. τ Q is the duration of the ramp that takes place in the time interval t ∈ [−τ Q /2, τ Q /2]. Using the terminology of Ref. [29], we distinguish protocols with λ i and λ f corresponding to different phases of the model (ramp crossing the critical point), and protocols with λ f = 0 (ramp ending at the critical point). We are going to refer to these two choices as transcritical protocol (TCP) and end-critical protocol (ECP), respectively. Certain observables exhibit markedly different behaviour depending on the protocol [38], hence both of them are of interest.
Ramps along the free fermion line (h = 0) have been studied extensively, especially in the spin chain. The time evolution of the free fermion modes with different momentum magnitudes decouple and only modes of opposite momenta {k, −k} are coupled by the evolution equation. One can make progress either by invoking the Landau-Zener description of transitions between energy levels or by numerically solving the set of two differential equations. Even analytical solutions are known for various ramp profiles [24,52]. These solutions can be simply generalised to the continuum field theory, providing us with an analytical tool to examine the KZ scaling and offering a benchmark for our numerical method. We refer the reader to Appendix B for the details.
The Kibble-Zurek mechanism is much less studied along the other integrable axis M = 0. As we noted above, in this direction ν = 8/15, so the KZ scaling is modified with respect to the well-investigated free fermion case. Although the model is integrable, the time evolution cannot be solved analytically, which highlights the importance of the nonperturbative numerical method that exploits the conformal symmetry of the critical model: the Truncated Conformal Space Approach (TCSA). Nevertheless, standard KZ arguments rely only on typical energy and distance scales of the model, consequently they should apply regardless of the presence of interactions. The scaling arguments can be supported by the analysis of the exactly known form factors of the model in the context of the adiabatic perturbation theory, to which we turn now.

Adiabatic Perturbation Theory
The adiabatic perturbation theory (APT) is a standard approach to study the response to a slow perturbation [86,25]. It was first used to describe the universal dynamics of extended quantum systems in the vicinity of a quantum critical point in Ref. [7]. Ever since the framework has become more elaborate by exploring the parallelism between APT and the Kibble-Zurek mechanism and generalizing the arguments to a wider variety of scaling quantities in different models [16,21,17,22,19,23,38]. In particular, it has already been applied with success in an integrable field theory, the sine-Gordon model [17]. In our current work we carry out an analogous reasoning to explore the implications of the APT statements in the E 8 Ising Field Theory. To this end, let us briefly sketch the basic concepts and assumptions underlying the framework of adiabatic perturbation theory as well as introduce some notations. Our discussion is based on the presentation of Ref. [22].
Assume that we want to solve the time-dependent Schrödinger equation: Using the basis of eigenstates of H(t) that are going to be called instantaneous eigenstates |n(t) , we can expand the time evolved state with coefficients α n (t): where the dynamical phase factor Θ n (t) = t t i E n (t ) dt is already included. The initial condition is that at t i the system is in its ground state |0(t i ) . Substituting this Ansatz into Eq. (2.16) yields a system of coupled differential equations for the coefficients α n (t). The resulting system of equations can be solved approximately for α n (t) using a few assumptions. First, the explicit time-dependence of the Hamiltonian is due to a time-dependent coupling constant λ that couples to some perturbing operator V so H(t) = H 0 + λ(t)V . Second, λ(t) is a monotonous function of time, hence one can perform a change of variables, and it changes slowly (that is the adiabatic assumption) such thatλ → 0. Then the resulting expression can be expanded in terms of powers ofλ. Assuming there is no Berry phase, the result up to leading order inλ is with λ = λ(t), and the dynamical phase with respect to the coupling is Θ n (λ) = λ E n (λ )/λ dλ .
Note that the phase factor exhibits rapid oscillations in the limitλ → 0. This can be exploited to identify the two possibly dominant contributions of integral Eq. (2.19) in this limit. First, a non-analytic part that comes from the saddle point of the phase factor at a complex value of coupling λ. It is exponentially suppressed with the inverse of the ratė λ. Second, there are contributions coming from the boundaries of the integration domain which can be obtained by integrating by parts and keeping terms to first order inλ yields the result This contribution can be viewed as a switch on/off effect as it is the consequence of a nonsmooth start or end of the ramp: it is nonzero if the first time derivative of the coupling has a discontinuity at the initial or final times. Ifλ i,f = 0 then one has to go to higher orders.
In general, a discontinuity in the ath derivative brings about the scaling α ∝ τ −a Q with the time parameter of the ramp τ Q [24]. We consider linear ramps (cf. Eq. (2.15)) so higher derivatives disappear and the small parameter of the perturbative expansion is 1/τ Q . We remark that Eq. (2.20) can be modified if the energy difference in the denominator vanishes at some time instant along the process, in that case the dependence of α onλ is subject to change (cf. Eq. (2.25) for low-momentum modes if the gap is closed).
The applicability of adiabatic perturbation theory, strictly speaking, requires that the overlap between the time-evolved state and the instantaneous ground state remains close to 1 [86]. This, however, imposes a constraint on the probability to be in an excited state rather than on the density of excitations. On the other hand, for quantum many-body systems in the thermodynamic limit the physical criterion for a perturbative treatment is to be in a low-density state [19]. Given that the Kibble-Zurek mechanism predicts that densities decay as a power law of the time parameter τ Q , in the limit τ Q → ∞ the above approximations are justified and we can use Eq. (2.19) to examine the Kibble-Zurek scaling. This reasoning predicts the correct scaling exponents in the transverse field Ising chain for various quantities [22,38]. Let us illustrate how they work in the case of the density of defects n ex after a linear ramp The states of the Ising chain participating in the dynamics are products of zero-momentum particle pair states with momentum k, hence the defect density can be expressed as 4 where α k = α k (λ f ) is the coefficient of a particle pair state |k, −k given by Eq. (2.19). To investigate the dependence on τ Q it is practical to introduce the rescaled variables (2.22) to remove the 1/τ Q dependence from the exponent of Eq. (2.19). The heart of the APT treatment of KZ scaling lies at the observation that the matrix element and energy difference appearing in the expression of α n take the following scaling forms: with the asymptotic behaviour F (x) ∝ x z and G(x) ∝ x −1/ν as x → ∞. These considerations yield that (2.26) Eq. (2.25) is analysed in the limit τ Q → ∞. In that case the limits of the integral over η are sent to ±∞ and one has to check whether the resulting integral converges or not. Substituting Eqs. (2.23) and (2.24) one can perform the integral in (2.26) in the limit η ζ ν i,f to determine the asymptotic behaviour The criterion for convergence then is 2z + 2/ν > 1, or, equivalently ν 1+zν < 2 [22]. In the opposite case the integral is divergent, indicating that to discard the contribution from high-energy modes in the limit τ Q → ∞ is not justified. The scaling brought about by all energy scales is quadratic τ −2 Q due to the discontinuity ofλ, cf. Eq. (2.20). Consequently, the case of equality ν 1+zν = 2 distinguishes between the Kibble-Zurek scaling determined by the exponent of τ Q in Eq. (2.25) and the quadratic scaling.

Application to the Ising Field Theory
These are the key themes of adiabatic perturbation theory as applied to model the Kibble-Zurek mechanism. Now we are going to show that these considerations can be generalised to the two integrable directions of the Ising Field Theory. In the case of the free field theory the generalisation of the arguments above is straightforward and it yields the same result as for the free fermion Ising chain. To apply the reasoning to the E 8 integrable model requires a bit of extra work. The complications are mainly technical, details are presented in Appendix A. Here we would like to highlight the key assumptions of the arguments only.
There are several major differences between the free fermion and the E 8 field theory: the spectrum of the latter exhibits eight stable stationary particles, moving particle states are built up by combining particles of various species. As a result, there are multiple kinds of many-particle states in contrast to the pair of a single particle species in the free field theory. Interactions between particles modify the simple p n = 2πn/L quantisation rule of momenta in finite volume L, leading to a nontrivial density of states in momentum space. Eigenstates of the theory are asymptotic scattering states labelled by the relativistic rapidity variable ϑ that is related to the energy and momentum of particle j as E j = m j cosh ϑ j and p j = m j sinh ϑ j .
To investigate the Kibble-Zurek scaling in this model we make several simplifying assumptions. First, we consider low-density states which is justified in the limit τ Q → ∞. Apart from being a necessary assumption to use the framework of adiabatic perturbation theory, it sets the ground for our second assumption: that is, we assume that the contribution from one-and two-particle states contribute dominantly to intensive quantities such as the defect and energy density. In contrast to the free fermion case, the time-evolved state in the E 8 model includes contributions from multiparticle states that do not factorize exactly to a product of particle pairs. On the other hand, the many-particle overlap functions still satisfy the pair factorisation up to a very good approximation given that the energy density of the non-equilibrium state is low [87,80] compared to the natural scale set by the final mass gap. Intuitively, the essence of this approximation is that due to large interparticle distance, the interactions between particles located far from each other can be neglected. Hence, the contribution of genuine multiparticle states is proportional to the probability of more than two particles located within a volume related to the correlation length. For a low-density state this probability is indeed tiny, hence the pair factorization is a good approximation. This assumption is also verified by previous works modeling the non-equilibrium dynamics of the Ising Field Theory that show that time evolution after sudden quenches is dominated by few-particle overlaps in the regime of low post-quench density [79,82,88].
Based on these assumptions, we can show that the arguments of APT generalise to an interacting field theory as well. Let us sketch the derivation for the excess heat density w that can be expressed as We evaluate this expression by calculating the α n coefficients as given by Eq. (2.19) in finite volume and then take the L → ∞ limit. Taking into account the finite volume expression of matrix elements in the E 8 model, we find that one-particle states contribute to the energy density with the right KZ exponent τ − ν ν+1 Q (for details see Appendix A.1). To the best of our knowledge, this is the first case when the KZ scaling of one-particle states is investigated in adiabatic perturbation theory.
The contribution of a two-particle state with species a and b is going to be denoted w ab and reads where ϑ ab is a function of ϑ determined by the constraint that the state has zero overall momentum. To take the thermodynamic limit one has to convert the summation to an integral over rapidities. The key observation to proceed is that the effects of the interactions are of O(1/L) and disappear in the limit L → ∞. Consequently, one can change the integration variable such that it goes over momentum instead of the rapidity. From then on, the derivation is identical to the free fermion case, although one has to check whether the scaling forms (2.23) and (2.24) apply for the dispersion and matrix elements of the E 8 theory as well. Observing that ϑ = arcsinh(p/m a ) = arcsinh [p/ (c|λ| ν )] with some constant c, one can see that the former is trivially satisfied with the right asymptotic F (x) ∝ x z . The latter equation regarding the scaling and the high-energy behaviour of the matrix element also holds in general, as one can verify in the E 8 model (see Appendix A). Hence, as long as the initial assumptions of low energy and approximate pair factorisation are valid, the adiabatic perturbation theory predicts KZ scaling of intensive quantities in the E 8 theory as well. Let us remark that the perturbative calculations indicate that the KZ scaling applies to each contribution coming from any one-particle state and two-particle branch separately. That is a nontrivial statement since the spectrum of the E 8 field theory is a result of a bootstrap procedure relying heavily on delicate details of the interaction, however, these details are overlooked by a first order perturbative calculation. Although we expect that the summed contribution of one-and two-particle states to the energy density satisfies the KZ scaling (in line with the generic reasoning of Sec. 2.1), the much stronger statement of APT concerning the scaling behaviour of separate branches does not necessarily hold true. We can draw an analogy with the form factor series expansion calculation of the central charge, where the result of the sum over multiparticle states is fixed by the c-theorem, while the separate terms vary greatly due to the details of the interaction [89].
We note that in the current case the ambiguity arises from taking the L → ∞ limit, since strictly speaking the adiabatic perturbation theory is sensible only if the ground state overlap remains close to 1, which is impossible for a finite density state in the thermodynamic limit. Previous calculations within the APT framework illustrate that this condition can be relaxed when calculating intensive quantities [19,38], demanding a low-density time-evolved state instead of one with almost unity overlap with the instantaneous ground state. Although this approach successfully captures qualitative features of the KZ scaling, the above considerations indicate that one has to be careful as to what extent to draw conclusions from it.

Truncated Conformal Space Approach
After introducing the perturbative approach to model the scaling laws of the Kibble-Zurek mechanism in the Ising Field Theory, let us now address a non-perturbative numerical method that can be used to verify the arguments above by explicitly simulating the dynamics. In the following we turn our focus to the Truncated Conformal Space Approach and discuss the underlying principles and its operation.
Numerical methods that are based on truncating the Hilbert space have a long history of capturing equilibrium properties of field theories (see [65] for a review). In particular, twodimensional field theoretical models that are defined by perturbing a conformal field theory or free theory by relevant operators are amenable to a very efficient numerical treatment, called the Truncated Conformal Space Approach (TCSA) [63,64]. The essential idea of the method is to compute the matrix elements of the perturbing operators in the basis of the unperturbed theory in finite volume where the spectrum is discrete. The resulting Hamiltonian matrix is then made finite dimensional by truncating the basis, hence the name of the method. Recently, it has been applied with success to model the non-equilibrium dynamics of different theories, in particular the Ising Field Theory [79,82,88,84]. We dedicate this section to briefly introduce the method and set up some notation along the course.
To model the Kibble-Zurek mechanism in the Ising Field Theory we define the theory in a finite volume L using periodic boundary conditions, so the space-time covers an infinite cylinder of circumference L. The basis states used by TCSA are the energy eigenstates of the c = 1/2 conformal field theory on the cylinder. The truncation keeps only a finite set of states that diagonalise the conformal Hamiltonian H 0 by discarding those with energy larger than a given cut-off E cut . The exact finite volume matrix elements of the primary fields σ and ε can be constructed on this basis by mapping the cylinder to the complex plane where conformal Ward identities can be utilised. Perturbing the CFT opens a mass gap ∆ that can be used to express the Hamiltonian matrix H in a dimensionless form for numerical calculations: where l = ∆L is the dimensionless volume parameter, ∆ φ is the scaling dimension of the field φ = σ, ε with ∆ σ = 1/8 and ∆ ε = 1. Hereκ is the dimensionless coupling constant that characterises the strength of the perturbation. The ramping protocol is thus realised in TCSA by tuningκ linearly in the dimensionless time ∆ i t, where ∆ i is the mass gap at the initial time instant. All quantities are measured in appropriate powers of ∆ i along the course of the ramp. Referring to the different physical content of the theories that result from the choice of σ or ε we use different notation for the mass gap in this work. The σ perturbation yields the E 8 spectrum with eight stable particles hence the notation for the mass gap in this case is m 1 , the mass of the lightest particle. The ε direction corresponds to a free fermion field theory with a single species so we simply denote ∆ as m the mass of the elementary excitation.
The success of TCSA to model the physical theory without an energy cut-off relies on its capability to suppress truncation errors as much as possible. Achieving higher and higher cut-offs is computationally demanding but the contribution of high energy states can be taken into account through a renormalisation group (RG) approach [90,91,73,92,77,93,94]. The RG analysis predicts a power-law dependence on the cut-off. Here we use a simpler extrapolation scheme using the powers predicted by the RG analysis which improves substantially the results obtained using relatively low cut-off energies. We express the recipe for extrapolation in terms of the conformal cut-off level N cut that is related to the energy cut-off as N cut = L/(2π)E cut . One can show that the results for some arbitrary quantity φ at infinite cut-off are related to TCSA data as where the α φ < β φ exponents are positive numbers depending on the scaling dimension of the perturbation, the operator in consideration and those appearing in their operator product expansion. Ellipses denote further subleading corrections that decay faster as N cut → ∞. The details of the extrapolation in various cases are detailed in Appendix C.
With this we have finished reviewing the basic concepts in the Kibble-Zurek mechanism and in the Ising Field Theory. We have introduced the two main methods that we use to study it: the numerical method of TCSA for simulating the dynamics and the scaling arguments in the context of APT that predicts that for the KZ scaling the presence of interactions in the E 8 theory makes no difference. We have outlined the following claims: the scaling behaviour observed on the transverse field Ising chain does not change in the continuum limit and that the only modification needed for the interacting E 8 model is to take into account the different scaling exponent ν. Before putting these claims to test by calculating the dynamics of one-point functions and observing the statistics of excess heat, we investigate the dynamics of energy eigenstates along the ramp in order to sketch an intuitive picture of how the Kibble-Zurek mechanism can be understood at the most fundamental level.

Work statistics and overlaps
We aim to study the evolution of the quantum state during the ramp, including the non-adiabatic regime, in detail. Using the TCSA method, we have access to microscopic data, which allows us to investigate the details of the dynamics. There are many possible quantities to consider: the correlation length, excitation densities, etc. In this section we adopt another, more microscopic perspective: we observe how instantaneous eigenstates get populated in the course of the ramp, how the adiabatic behaviour breaks down and how excitations are created. Looking at the fundamental components that conspire to create the well-known KZ scaling in a wide variety of quantities provides us with an intuitive and visual picture about what happens during the regime when adiabaticity is lost.
To this end, we first solve the time-dependent Schrödinger equation: with the initial state |Ψ 0 chosen to be the ground state of the initial Hamiltonian H(−τ Q /2). Since momentum is conserved all along the ramp and the initial state is a zero-momentum state, |Ψ(t) is also a P = 0 state for all t.
To characterise how the energy eigenstates get populated we can generalise the statistics of work function [95] to each time instance along the course of the ramp, defining an instantaneous statistics of work function where the sum is running over the spectrum of the instantaneous Hamiltonian H(t) with eigenvalues E n (t) and eigenstates |n(t) . Here g n (t) are the overlaps of the time-evolved state with the instantaneous eigenstates: W is called to the total work performed by the non-equilibrium protocol. P (W , t) is nonzero only ifW ≥ E 0 (t) − E 0 (0). In the following we focus only on the statistics of the In order to draw a clear picture of what happens for ramps within the reach of KZM, we present the two sections of P (W, t): first, only the |g n (t)| 2 overlap amplitudes with respect to time and second, the snapshot of P (W, t) at some time instant t.

Ramps along the free fermion line
Let us start with the exactly solvable dynamics, i.e. the free fermion line of the model (2.11) corresponding to h = 0. The time-dependent coupling is the free fermion mass, λ(t) = M (t). Our ramp protocol is a simple linear ramp profile that is symmetric around the critical point: where M i is the initial value of the coupling at t = −τ Q /2. As discussed in Sec. 2.2, the critical exponents in this case are ν = 1, z = 1, so the Kibble-Zurek time (2.2) scales as τ KZ ∼ √ τ Q . For testing the various scaling forms we need to have a specified value of τ KZ which we simply set as where m = |M i | is the mass gap at the start of the ramp. Depending on the sign of M i , the ramp is either towards the ferromagnetic phase or the paramagnetic phase; we are going to present our results in this order.

The paramagnetic-ferromagnetic (PF) direction
Ramps starting from the paramagnetic phase are defined by M i > 0. In this case the ground state is non-degenerate and lies in the Neveu-Schwarz sector, so the time evolved state is orthogonal to the Ramond sector subspace for all times (see Sec. 2.2).
Analogously to the lattice dynamics, starting from the ground state at a given M i , only states consisting of zero-momentum particle pairs have nonzero overlap with the time evolved state, moreover, the different pairs of momentum modes {p, −p} decouple completely. In finite volume L the momentum is quantised as p n = 2πn/L, where n is half-integer in the NS sector. To solve the dynamics we follow the approach of [52] and use the Ansatz: where |0 p,t and |1 p,t denote the instantaneous ground and excited states of the two-level system at time t along the ramp. The coefficients a p (t) and b p (t) satisfy |a p (t)| 2 +|b p (t)| 2 = 1 and they can be expressed via the solutions of two coupled first order differential equations (for details see the Appendix B). The population of mode p is given by n p (t) = |b p (t)| 2 .
Although the equations can be solved exactly, numerical integration is more suitable for our purposes. Hence, strictly speaking, referring to this solution as 'analytical' is not entirely precise. From now on, when we use the term 'analytical' we mean the 'numerically exact' procedure outlined above.
Apart from this solution of the dynamics, we can calculate the population of energy eigenstates numerically with TCSA. This is a benchmark for our numerical method as it is contrasted with a numerically exact calculation. We can compare Eq. (3.6) with Eq. (3.3) to express the overlap g of a state which consists of only a single particle pair with momentum p: where the product goes over the infinite set of quantised momenta in finite volume. It is straightforward to generalise Eq. (3.7) to express the overlap of any state with the pair structure of the free spectrum with the time-evolved state.
In practice, we truncate this product at some finite p max , since the goal is to match the analytic results with TCSA that operates with a truncation of its own. The one-mode cut-off of the analytic method and the many-body cut-off of TCSA cannot be brought to one-to-one correspondence with each other. However, overlaps are very sensitive to the number of states kept in each expansion, due to the constraint n |g n | 2 = 1. Hence, our choice for the energy cutoff of TCSA for these figures is motivated by the goal to have the best possible match of the two approaches. Note that this is a single parameter for all the states.
The time evolution of the overlaps is presented in Fig. 3.1. Dots correspond to the solution of the differential equations for each mode and continuous lines denote TCSA data obtained by solving the many-body dynamics numerically. Fig. 3.1a depicts a curious behaviour of the second largest overlap in TCSA: the corresponding line seemingly consists of many different segments. This is a consequence of level crossings and the errors of numerical diagonalisation near these crossings. The state in question consists of two twoparticle pairs and as the mass scale M is ramped its energy increases steeper than that of high-momentum states with only a single pair, hence the level crossings. At each crossing the numerical diagonalisation cannot resolve precisely levels in the degenerate subspace, so the resulting overlap is not accurate. This accounts for the most prominent difference between the numerical and analytical results. Apart from that, the agreement is quite satisfactory.
The light green background corresponds to the naive impulse regime t ∈ [−τ KZ , τ KZ ]. Of course this is only a crude estimate for the time when adiabaticity breaks down as Eq. (3.5) is strictly valid only as a scaling relation. Nevertheless, most of the change in each state population indeed happens within this coloured region. This statement is even more accentuated by Fig. 3.1b, that is, for a slower ramp. Comparing the two panels of Fig. 3.1 we observe that increasing the ramp time the probability of adiabaticity increases while  We can also calculate the energy resolved version of the above figures, i.e. the instantaneous statistics of work, P (W, t). We present this quantity in Fig. 3.2. The different ridges correspond to "bands" of 2-particle, 4-particle etc. states with energy thresholds E = 2M, 4M, . . . . The ridges diverge linearly in time, displaying the linear dependence of the gap on the linearly tuned M coupling. This figure illustrates the validity of the KZ arguments: low-energy bands dominate the excitations, and in each band, the modes with the lowest momenta (longest wavelengths) near the thresholds are the most prominent. This feature is similar to what was observed on the lattice in Ref. [37].

The ferromagnetic-paramagnetic (FP) direction
The ferromagnetic ground state is twofold degenerate in infinite volume. For the initial state we choose the state with maximal magnetisation corresponding to the infinite volume symmetry breaking state: |Ψ 0 = 1 √ 2 (|0 R + |0 NS ). As both sectors are present in the initial state, the time-evolved state also overlaps with both sectors. This provides yet another benchmark for our numerical approach and also a somewhat richer landscape of the overlap functions.
As one can see in Fig. 3.3, the dynamics are very similar to the PF case with the main difference coming from the fact that both sectors contribute. The different behaviour of the two vacua stems from the different available momentum modes in each sector: in the Ramond sector the momenta are larger in the lowest available modes and consequently   they are less likely to be excited.

Ramps along the E 8 line
After investigating the free fermion line, we now turn to the behaviour of overlaps in the other integrable direction, i.e. for ramps along the E 8 axis defined by the protocol   for t ∈ [−τ Q /2, τ Q /2]. The scaling dimension of the perturbing operator σ is ∆ σ = 1/8, so critical exponent ν is different in this direction from the free fermion case: ν = 1/(2−∆ σ ) = 8/15 (cf. Eq. (2.14)). This implies that the Kibble-Zurek time (2.2) is given by where, similarly to the free fermion case, the choice of the proportionality factor being 1 is just a convention. Let us first take an overview of the dynamics by looking at the time-dependent work statistics P (W, t) shown in Fig. 3.4. Notice that in accordance with the Kibble-Zurek scenario, predominantly low-energy and low-particle modes get excited in the course of the ramp. In the E 8 theory with multiple stable particles, the time evolved state has finite overlap not only with states consisting of pairs but also with states containing standing particles with zero momentum, including multiparticle states with a single such particle. We can observe that the energy distribution has peaks at some finite energy values, but low-momentum modes dominate for all branches (denoted by dashed lines of the same colour). This can be seen more clearly in Fig. 3.5 which presents P (W ) at the end of two ramps that differ in duration. Solid vertical lines indicate the energies of states consisting of standing particles only, i.e. combinations of particle masses.
As discussed at the end of Sec. 2.3 (and derived in detail in App. A), perturbation theory predicts that the overlaps of these standing particle states decay uniformly with the quench time as τ −8/23 Q . Fig. 3.5 clearly illustrates that this is not the case: as the average excess heat diminishes, the overlap of low-lying states increase instead of decrease. However, as we are going to show later, both quench times are within the KZ scaling region and the scaling of the excess heat does satisfy Eq. (2.6).

Probability of adiabaticity
To study the Kibble-Zurek scaling using the TCSA, it is important to identify the time scale on which it is valid. For a finite volume method the time scale is limited from above by the onset of adiabaticity (cf. Eq. (2.9)) and also from below due to the natural time scale of the theory that is related to the mass gap before and after the ramp. A control quantity that can be used to fix the domain of τ Q where the Kibble-Zurek scaling applies is the probability to be adiabatic after the ramp, P (0, t f ). This overlap is exponentially suppressed with the volume, but its logarithm is proportional to the density of quasiparticles n ex , such that − log(P (0))/L ∝ n ex . Within the domain of validity for the Kibble-Zurek scaling the density scales according to Eq. (2.4), i.e. decays as a power law with τ Q . However, at the onset of adiabaticity it is exponentially suppressed [6,13]. To explore the time scale mentioned above connected to volume parameters available for our calculation, we investigate the logarithm of the ground state overlap P (0) after the ramp.
For ramps along the free fermion line there are two ways to evaluate P (0). The first follows from the numerically exact solution of the problem in the scaling limit (see Appendix B). Second, we can use TCSA to calculate the ground state overlap. The onset of adiabaticity occurs at different quench times τ Q depending on the volume parameter. Then the claim that for a given volume L we can observe the KZ scaling -as opposed to adiabatic behaviour -can be supported by the observation that changing the volume does not alter the KZ scaling. Fig. 3.6a presents the comparison of the two methods with the slope of the KZ scaling as a guide to the eye. Apart from the very fast ramps, the two methods coincide with each other. We note that the onset of adiabaticity signalled by the strong deviation of different volume curves from each other and from the τ −1/2 Q line is not an abrupt change but rather a smooth crossover. Nevertheless, we can identify that/ for mτ Q ≈ 5 · 10 0 . . . 10 2 the Kibble-Zurek scaling is satisfied to a good precision using the volume parameters available to the numerical method.
In the E 8 model we can only resort to the results of TCSA. Fig. 3.6b shows that the logarithm of the ground state overlap scales as the density of quasiparticles for large enough τ Q . Although the KZ scaling sets in later, i.e. for larger τ Q than in the free fermion case, it is persistent up to the maximum ramp duration available to our numerical method. This is due to the fact that the exponent appearing in Eq. (2.9) is larger for the E 8 model and consequently the onset of adiabaticity occurs for a slower ramp in the same volume.

Ramps ending at the critical point
As detailed in Section 2.3, we expect the generic scaling arguments of APT for the Kibble-Zurek mechanism (Eqs. (2.23) and (2.24)) to be valid for ramps along both inte-  grable lines of the model. A direct consequence of this claim is that the high-energy tail of the function |K(η)| 2 decays as η β with β = −2z − 2/ν (cf. Eq. (2.27)). This behaviour is important in view of the convergence properties of the integrals of the form (2.25).
To investigate the decay of high-energy overlaps with TCSA, we consider ramp protocols along the two integrable lines of the parameter space that end at the conformal point (ECP ramps). There are two reasons for this choice of protocol: first, TCSA uses the conformal basis and hence expected to be the most accurate at the critical point. Second, the dispersion relation is E(k) = |k| in this case, so the high-energy tail of P (W ) decays with the same power law as |α(k)| 2 . Since k and η are related by a simple rescaling with the appropriate power of τ Q , the high energy tail of P (W ) should decay as W β at the critical point as far as the perturbative approach is correct, i.e. for slow enough ramps.
On the free fermion line we have z = ν = 1, so β = −2z − 2/ν = −4, while for an E 8 ramp ν = 8/15 and the predicted exponent of the decay is β = −23/4. We remark that this can be contrasted with the high-energy tail of pair overlaps for sudden quenches. For quenches along the free fermion line the exact solution yields β = −2 [96,97,79], while in the E 8 model the high energy tail of the perturbative expression decays with β = −15/4 [88], so β = −2/ν in both cases. The additional term of −2z is the result of the adiabatic driving which suppresses the excitation of high energy modes.
In Fig. 3.7 we present the TCSA data and the slope of the straight line fitted to the logarithmic data. The two exponents are well separated and captured approximately correctly by the data. We note that Fig. 3.7a is analogous to Fig. 2c of Ref. [37] that reported a W −8 decay. This is at odds with the prediction deduced from generic scaling arguments using APT and also with our TCSA results that favor the β = −4 exponent. Fig. 3.7 is in agreement with the numerous observations [7,16,24,38] that adiabatic perturbation theory captures the correct Kibble-Zurek scaling in the free fermion theory and demonstrates that it applies also in the interacting E 8 integrable model. This is evidence that the arguments of APT can be generalised to this nontrivial theory which in turn implies that the Kibble-Zurek scaling can be observed there as well.

Dynamical scaling in the non-adiabatic regime
In this section we explore the dynamical scaling aspect of the Kibble-Zurek mechanism in the Ising Field Theory considering two one-point functions. We focus on the energy density and the magnetisation, both of which are important observables in the theory.
The energy density over the instantaneous vacuum or the excess heat density is defined as where the Hamiltonian H(t) has an explicit time dependence governed by the ramping protocol and E 0 (t) is the ground state of the instantaneous Hamiltonian H(t). In accordance with Eq. (2.6), the excess heat for different ramp rates is expected to collapse to a single scaling function: where d = 1 is the spatial dimension, ∆ H = z is the scaling dimension of the energy and the second equation follows from τ KZ = ξ z KZ . For ramps along the free fermion line the energy density can be obtained from the solution of the exact differential equations using the mapping to free fermions, yielding essentially exact results.
The magnetisation operator σ that corresponds to the order parameter has scaling dimension ∆ σ = 1/8 hence is expected to satisfy the following scaling in the impulse regime (z = 1): In contrast to the energy density, the magnetisation is much harder to calculate even in free fermion case as it is a highly non-local operator in terms of the fermions.

Free fermion line
We start with the free fermion line where exact analytical results are available. In Fig.  4.1a we observe the scaling behaviour (4.2) for several ramps from the paramagnetic to the ferromagnetic phase. Both the analytic calculations and the TCSA data, extrapolated in the cutoff, retain the scaling and the numerics agree almost perfectly with the exact results. The inset shows that the non-rescaled curves deviate substantially from each other.
As Fig. 4.1a shows, the collapse of the curves is perfect even well beyond the end of the non-adiabatic regime, in agreement with the observation and arguments of Ref. [31]. This can be understood in view of the eigenstate dynamics presented in Sec. 3. The relative population of energy eigenstates does not change substantially in the post-impulse regime and the increase in energy density then is merely due to the increasing gap ∆(t) as the coupling is ramped. The energy scale increases identically for all quench rates which in turn leads to the collapse of different curves. This argument can be formalised for the general setup of Sec. 2.1 as where n ex is the density of defects that is constant well beyond the impulse regime and scales as τ −d/z KZ . The gap scales as (t/τ Q ) zν and we used that (τ KZ /τ Q ) aνz ∝ τ −1 KZ . The result shows that w(t τ KZ ) is a function of t/τ KZ . In the present case a = ν = z = 1, which explains the linear behaviour seen in Fig. 4.1a.
The scaling behaviour of the magnetisation (4.3) is checked in Fig. 4.1b. The scaling is present most notably in terms of the frequency of the oscillations beyond the non-adiabatic window. Due to truncation errors of the TCSA method (see Appendix C), the predicted scaling is not reproduced perfectly in terms of the amplitudes and neither in the first half of the non-adiabatic regime. This is also the reason why the various curves do not collapse perfectly for times t < −τ KZ where the scaling should also hold according to Eq. (2.7).
The frequency of the late time oscillations is increasing with time. The oscillations can  be fitted with the function f (t) = A cos [m(t) · t + φ] which demonstrates that the oscillations originate from one-particle states whose masses and thus the frequency increases in time with the gap. We remark that this is analogous to sudden quenches in the Ising Field Theory where the presence of one-particle oscillations is supported by analytical and numerical evidence [96,79,82]. The oscillations appear undamped well after the impulse regime t/τ KZ 1. We remark that for sudden quenches the decay rate of the oscillations depends on the post-quench energy density [97,96]. We expect the same to apply for ramps as well, but here the energy density is suppressed for slower ramps so the damping cannot be observed during a finite ramp. In contrast, the decay of oscillations in the dynamics of the order parameter after the ramp is observed in Ref. [37] in the spin chain.

Ramps along the E 8 axis
The dynamical scaling is well understood for the free fermion model on the lattice, and in the previous sections we demonstrated that they apply in the continuum scaling limit as well. The same aspect of the other integrable direction of the Ising Field Theory is yet unexplored. We now present how the simple scaling arguments of the KZM apply in a strongly interacting model. The dynamics in the E 8 model cannot be treated exactly due to the interactions but the numerical method of TCSA can be applied to simulate the time evolution. Truncation errors are expected to be less substantial since the σ perturbation of the CFT is more relevant and exhibits faster convergence compared to the free fermion model (cf. Fig. 3.7). Hence using the conformal eigenstates as a basis of the Hilbert space is expected to be a better approximation.
As discussed above, the scaling is modified compared to the free fermion model due to the different exponent ν = 8/15, so the Kibble-Zurek time scale τ KZ depends on the ramp time τ Q as τ KZ = τ 8/23 Q . We demonstrate this scaling in the following for the dynamics of the energy density and the magnetisation.
Let us first discuss the scaling of the energy density presented in Fig. 4.2a. Similarly to the free fermion case, one observes an almost perfect collapse of the curves after crossing the critical point, and the collapse is sustained beyond the impulse regime where now Eq. Note that the above argument relies on the fact that the scaling properties of the energy density can be determined by considering it as the product of some defect density and a typical energy scale. For more complex quantities, such as the magnetisation for example, a similar argument does not apply, as Fig. 4.2b demonstrates. The curves deviate after the non-adiabatic regime but the collapse in the early adiabatic regime is perfect.

Cumulants of work
So far we have gained insight in the KZM by examining the instantaneous spectrum directly and demonstrated the relevance of the Kibble-Zurek time scale in dynamical scaling functions of local observables. In this section we aim to demonstrate that the Kibble-Zurek scaling is present in an even wider variety of quantities: the full statistics of the excess heat (or work) during the ramp is subject to scaling laws of the KZ type as well.
A particularly interesting result of the free fermion chain (already tested experimentally, cf. Ref. [46]) is that apart from the average density of defects and excess heat, their full counting statistics is also universal in the KZ sense: all higher cumulants of the respective distribution functions scale according to the Kibble-Zurek laws [34,38]. The scaling exponents depend on the protocol in the sense that they are different for ramps ending at the critical point (ECP) and those crossing it (TCP). As Ref. [35] demonstrates, the universal scaling of cumulants can be observed in models apart from the transverse field Ising spin chain, hence it is natural to explore their behaviour in the Ising Field Theory.
The cumulants of excess work are defined via a generating function ln G(s): where the expectation value is taken with respect to the time-evolved state. The cumulants κ i are the coefficients appearing in the expansion of the logarithm: The first three cumulants coincide with the mean, the second and the third central moments, respectively. Assuming that the generating functions satisfy a large deviation principle [99,38], all of the cumulants are extensive ∝ L. Consequently, we are going to focus on the κ i /L cumulant densities. Elaborating on the framework of adiabatic perturbation theory presented in Sec. 2.3, we can argue that the scaling behaviour of the cumulants of the excess heat are not sensitive to the presence of interactions in the E 8 model and take a route analogous to Ref. [38] to obtain the KZ exponents. The core of the argument is the following: the Kibble-Zurek scaling within the context of APT stems from the rescaling of variables (2.22) which yields Eq. (2.25) from Eq. (2.21). The rescaling concerns the momentum variable that originates from the summation over pair states. Now consider that cumulants can be expressed as a polynomial of the moments of the distribution: where λ = {n 1 , n 2 , . . . , n k } is a partition of the integer index n with |λ| = k ≥ 2, and α λ are integer coefficients. The moments are defined for the excess heat as Let us note that the integration variable subject to rescaling in Eq. (2.22) originates from taking the expectation value. Consequently, in the limit τ Q → ∞ terms consisting of powers of lower moments are suppressed compared to µ n , because they are the product of multiple integrals of the form (2.25). So the scaling behaviour of κ n equals that of µ n , which is defined with a single expectation value, hence its scaling behaviour is given by the calculation in Sec. 2.3. We remark that this line of thought is completely analogous to the arguments of Ref. [38]. According to the above reasoning, all cumulants of the work and quasiparticle distributions in the E 8 model should decay with the same power law as τ Q → ∞.
To put the claims above to test, we follow the presentation of Ref. [38] and we discuss the two different scaling for the cumulants: first considering ramps that end at the critical point then examining ramps that navigate through the phase transition.

ECP protocol: ramps ending at the critical point
For ramps that end at the critical point one may apply the scaling form in (2.6) since the final time of such protocols corresponds to a fixed t/τ KZ = 0. The resulting naive scaling dimension of a work cumulant κ n is then easily obtained since it contains the product of n Hamiltonians with dimension ∆ H = z = 1. Consequently, we expect where we used Eq. (2.2). However, the arguments of adiabatic perturbation theory [38] as outlined in Sec. 2.3 demonstrate that this naive scaling is true only if the corresponding quantity is not sensitive to the high-energy modes. However, using APT one can express the cumulants similarly to the defect density in Eq. (2.25). If the corresponding rescaled integral does not converge that means the contribution from high-energy modes cannot be discarded and the resulting scaling is quadratic with respect to the ramp velocity: τ −2 Q . The crossover happens when aν(d + nz)/(aνz + 1) = 2; for smaller n the KZ scaling applies while for larger n quadratic scaling applies with logarithmic corrections at equality [22].
For the free fermion line ν = 1 (a = d = z = 1) and the crossover cumulant index is n = 3. Fig. 5.1 justifies the above expectations for the three lowest cumulants by comparing the numerically exact solutions to TCSA results. TCSA is most precise for moderately slow quenches and the first two cumulants. There is notable deviation from the exact results in the case of the third cumulant although the scaling behaviour is intact. The deviation does not come as a surprise since the fact that the integral of adiabatic perturbation theory does not converge means that there is substantial contribution from all energy scales including those that fall victim to the truncation. Fig. 5.1 also demonstrates that for very slow quenches finite size effects can spoil the agreement between exact results and TCSA. This is the result of the onset of adiabaticity (cf. Fig. 3.6a).
We expect identical scaling behaviour from the other integrable direction of the Ising Field Theory in terms of τ KZ that translates to a different power-law dependence on τ Q . Indeed this is what we observe in Fig. 5.2. In this case there is no exact solution available hence solid lines denote the expected scaling law instead of the analytic result. The figure is indicative of the correct scaling although finite volume effects are more pronounced as the duration of the ramps is larger than earlier.

TCP protocol: ramps crossing the critical point
For slow enough ramps that cross the critical point and terminate at a given finite value of the coupling which lies far from the non-adiabatic regime where (2.6) applies, the excess work density scales identically to the defect density. This is due to the fact that the gap that defines the typical energy of the defects is the same for ramps with different τ Q and the excess energy equals energy scale times defect density. It is demonstrated in Ref. [38] that higher cumulants of the excess work share a similar property: their scaling dimension coincides with that of the mean excess work, consequently all cumulants of the defect number and the excess work scale with the same exponent. As we argued above, this claim is expected to be more general than free theories and in particular we claimed that it holds in the E 8 model. Fig. 5.3 demonstrates the validity of this statement for the second cumulant. In line with the reasoning presented earlier (cf. Eq. (5.3) and below), the subleading terms are more prominent than in the case of the first cumulant (the excess heat) and KZ scaling is observable only for larger τ Q . Higher cumulants do not exhibit the same scaling within the quench time window available for TCSA calculations. Due to the increasing number of terms in the expressions with moments for the nth cumulant κ n , we expect that the Kibble-Zurek scaling occurs for larger and larger τ Q , on time scales that are not amenable to effective numerical treatment as of now. Nevertheless, the behaviour of the second cumulant still serves as a nontrivial check of the assumptions that were used in Sec. 2.3 to apply APT to the E 8 model. As the argumentation did not rely explicitly on the details of the interactions in the E 8 theory, rather on the more general scaling behaviour of the gap (2.23) and the matrix element (2.24), we expect that a similar behaviour of the cumulants is observable in other interacting models exhibiting a phase transition.

Conclusions
In this paper we investigated the Kibble-Zurek scaling in the context of continuous quantum phase transitions in the Ising Field Theory. The KZ scaling describes the universal dependence of a range of observables on the quench rate and it is connected to the breakdown of adiabatic behaviour due to a critical slowing down near the phase transition. The Ising Field Theory accommodates two types of universality in terms of the static critical exponent ν that corresponds to two integrable models for a specific choice of parameters in the space of couplings. One of them describes a free massive Majorana fermion and it exhibits a completely analogous KZ scaling to the transverse field Ising chain that can be mapped to free fermions. Building on the lattice results, we expressed the nonequilibrium dynamics through the solution of a two-level problem and explored the Kibble-Zurek mechanism in terms of instantaneous eigenstates and various observables, including local operators and cumulants of the distribution of the statistics of work.
We have shown that the adiabatic-impulse-adiabatic scenario is qualitatively correct at the most fundamental level of quantum state dynamics. That is, in the sense that we can identify a non-adiabatic "impulse" regime where the most substantial change in the population of eigenstates happens, preceded and followed by a regime of adiabatic dynamics where these populations are approximately constant. We demonstrated that the relative length of the impulse regime compared to the duration of the ramp decreases as the time parameter of the ramp τ Q increases. This decrease happens according to the scaling forms dictated by the Kibble-Zurek mechanism. Although this simple picture has been put to test from many aspects in earlier works, the observation that it applies at the fundamental level of quantum states is still noteworthy.
We established parallelisms between the lattice and continuum dynamics for an extended set of scaling phenomena from the dynamical scaling of local observables to the universal behaviour of higher cumulants of the work. These analogies do not come as a surprise but their analysis in a field theoretical context is a novel result. Apart from generalizing recently understood phenomena on the lattice to the continuum, these observations serve as a benchmark for our numerical method, the Truncated Conformal Space Approach. Comparing with analytical solutions available in the free fermion theory, we have illustrated the capacity of this method to capture the intricate quantum dynamics behind the Kibble-Zurek scaling near quantum critical points. In spite of operating in finite volume, it is capable of demonstrating the presence of scaling laws within a wide interval of the time parameter τ Q without substantial finite size effects. This is of paramount importance in the demonstration that the KZ scaling is not limited to the noninteracting dynamics within the Ising Field Theory.
The second integrable direction in the coupling space of the IFT corresponds to the famous E 8 model with its affluent energy spectrum exhibiting eight stable particle states. One of the essential results of our work is that the Kibble-Zurek mechanism is able to account for the universal scaling of this strongly interacting model near the quantum critical point. In order to have a solid case for this observation, we elaborated on the framework of adiabatic perturbation theory and applied its basic concepts to the E 8 model. While a refined version of the originally suggested adiabatic-impulse-adiabatic scenario predicts universal dynamical scaling of local observables in the non-adiabatic regime (which we also verified using TCSA, see Sec. 4), employing APT to address the nonequilibrium dynamics provides perturbative arguments for the universal scaling of the full counting statistics of the excess heat and number of quasiparticles. This reasoning has been used recently to explain the universal scaling of work cumulants in a free model [38]. In this work we have taken the next step and discussed its implications for the interacting E 8 field theory. We argued that the interactions do not alter the universal scaling of cumulants and demonstrated this in Sec. 5 for the first cumulants both for end-critical and transcritical ramp protocols. We remark that our argument is in fact quite general and mostly relies on the small density induced by the nonequilibrium protocol. Since the KZ scaling predicts that the dynamics is close to adiabatic as τ Q → ∞, this is a sensible assumption. Consequently, the result is expected to hold generally, i.e. all cumulants of the excess work should scale with the scaling exponents predicted by adiabatic perturbation theory irrespective of the interactions in the model.
We note that there are several possible future directions. It is particularly interesting to test the scaling behaviour of "fast but smooth" ramps versus sudden quenches in the coupling space of field theoretical models [100,101,102]. The presence of universal scaling at fast quench rates is remarkable though to implement an infinitely smooth ramp in an interacting theory that is not amenable to exact analytic treatment is not trivial. Another fruitful direction to take is the exploration of nonintegrable regimes within the Ising Field Theory and examine the interplay between the physics related to integrability breaking and the Kibble-Zurek scenario. Our findings suggest that the latter is in fact quite general but its validity in a generic non-integrable scenario remains to be tested.

A Application of the adiabatic perturbation theory to the E 8 model
To use the framework of adiabatic perturbation theory in the E 8 model we assume that the time-evolved state can be expressed as with the dynamical phase factor Θ n (t) = t t i E n (t ) dt . We also assume that there is no Berry phase and thus to leading order in the small parameterλ the α n coefficients take the form Higher derivatives as well as higher order terms inλ are neglected from now on. The α n coefficients can be used to formally express quantities that have known matrix elements on the instantaneous basis of the Hamiltonian: In what follows, we present the evaluation of this sum -approximately, under conditions of low energy density discussed in the main text -for the case of O(t) = H(t) − E 0 (t) in the E 8 model. To generalise this calculation to the defect density or to higher moments of the statistics of work function is straightforward. The work density (or excess heat density) after the ramp reads The spectrum of the model consists of 8 particle species A a , a = 1, . . . , 8 with masses m a . The energy and momentum eigenstates are the asymptotic states of the model labelled by a set of relativistic rapidities {ϑ 1 , ϑ 2 , . . . ϑ N } and particle species indices {a 1 , a 2 , . . . a N }: with energy E n = N i=1 m a i cosh(ϑ i ) and momentum p n = N i=1 m a i sinh(ϑ i ). The summation in Eq. (A.4) in principle goes over the infinite set of asymptotic states. As discussed in the main text, for low enough density we can approximate the sum in Eq. (A.4) with the contribution of one-and two-particle states, analogously to the calculation in the sine-Gordon model in Ref. [17].
A.1 One-particle states Contribution of the one-particle states can be expressed as where m a is the mass of the particle species a and the summation runs over the eight species. We can write the coefficient α a as where {0} a (λ)| denotes the asymptotic state with a single zero-momentum particle. The matrix elements and masses depend on λ through the Hamiltonian that defines the spectrum. The matrix element can be evaluated as For an E 8 ramp that conserves momentum, V is the integral of the local magnetisation operator σ(x): V = L 0 σ(x)dx. Utilizing this we further expand where the square root in the denominator emerges from the finite volume matrix element [103] and F σ a is the (infinite volume) one-particle form factor of the magnetisation operator. It only depends on the coupling λ through its proportionality to the vacuum expectation value of σ. The particle masses scale as the gap: m a (λ) = C a |λ| zν , where C a are some constants. This allows us to write (A.10) We can perform the integral in the exponent that leads to a τ Q |λ| 1+zν dependence there.
To get rid of the large τ Q factor in the denominator, we introduce the rescaled coupling ζ with ζ = λτ whereC a and C a are constants that depend on C a , the one-particle form factors and the critical exponents. We note the integral is convergent for large ζ due to the strongly oscillating phase factor and also for ζ → 0 since 2ν − 1 − 3/2zν = −11/15 in the E 8 model. Substituting z = 1 in the exponent of τ Q leads to the correct KZ exponent of a relativistic model, ν/(1 + ν).
A.2 Two-particle states The contribution of a two-particle state with species a and b is going to be denoted w ab and reads where ϑ ab is a function of ϑ determined by the constraint that the state has zero overall momentum. The summation goes over the rapidities that are quantised in finite volume L by the Bethe-Yang equations: where I i are integers numbers and is the scattering phase shift of particles of type a and b. For a two-particle state Eq. (A.14) amounts to two equations of which only one is independent due to the zero-momentum constraint. It readsQ (ϑ) = m a L sinh ϑ + δ ab (ϑ − ϑ ab ) = 2πI , I ∈ Z . (A. 16) In the thermodynamic limit L → ∞ the summation is converted to an integral with the integral measure dϑ 2πρ (ϑ), whereρ(ϑ) is the density of zero-momentum states defined bỹ where Φ(ϑ) is the derivative of the phase shift function. The resulting integral is The α ϑ (λ f ) term can be expressed as (cf. Eq. (A.2) (A.19) Analogously to the one-particle case we can evaluate the matrix element in the E 8 field theory as where F σ ab (ϑ 1 , ϑ 2 ) is the two-particle form factor of operator σ in the E 8 field theory and the density factor is the Jacobian of the two-particle Bethe-Yang equations (A.14) arising from the normalisation of the finite-volume matrix element [103]. It can be expressed as Observing Eqs. (A.17) and (A.21) one finds that the details of the interaction enter via the derivative of the phase shift function but crucially, they are of order 1/L compared to the free field theory part. So leading order in L we find that A change of variables in the outer integral to the one-particle momentum p = m a sinh ϑ we obtain Now we can introduce the momentum p in the inner integral as well by noting that the energy can be expressed as a function of momentum via the relativistic dispersion and that the relativistic rapidity also ϑ = arcsinh(p/m). Since m ∝ |λ| zν with z = 1 any expression that is a function of ϑ can be expressed as a function of p/|λ| ν . Having this in mind, the result is analogous to the free case so all the machinery developed there can be used. The key assumptions from this point regard the scaling properties of the energy gap and the matrix element G(ϑ) in this brief notation: These equations are trivially satisfied with the proper asymptotics for F (x) ∝ x z . For G(x) one can verify using that in the E 8 model we have where we neglected the O(1/L) term from the finite volume normalisation and used σ ∝ λ 1/15 , m ∝ λ 8/15 . F ab (ϑ, ϑ ab ) is the two-particle form factor of the E 8 theory that does not depend on the coupling. They satisfy the asymptotic bound [89]: Since the matrix elements considered here are of zero-momentum states, ϑ → ∞ means ϑ ab → −∞ and F σ ab (ϑ, ϑ ab ) ≤ exp(∆ σ ϑ) as the form factors depend on the rapidity difference. Dividing by the factor exp(2ϑ) in the denominator yields the correct asymptotics G(x) ∝ x ∆−2 = x −1/ν as an upper bound due to Eq. (A.27). We remark that the scaling forms (A.24) hold true for any value of the coupling λ in the field theory, in contrast to the lattice where they are valid only in the vicinity of the critical point. From this perspective Eq. (A.24) follows from the definition of the field theory as a low-energy effective description of the lattice model near its critical point.
As a consequence, one can introduce new variables in place of λ and p such that the explicit τ Q dependence disappears from the integrand. This is achieved by the following rescaling: (A.28) The result for the energy density is (A. 29) In terms of scaling there are two options: first, let |λ f | = 0 hence ζ f → ∞ in the KZ scaling limit τ Q → ∞. Then the energy gap at p → 0 is a constant and E p=0 (λ f ) can be brought in front of the integral. If it converges, Eq. (A.29) completely accounts for the KZ scaling.
Second, if |λ f | = 0, the energy gap is E p ∝ p z and an additional factor of τ − ν 1+zν Q appears in front of the integral. Note that this is the scaling of κ 1 on Fig. 5.2. The high-energy tail of the integrand is modified due to the extra term of η z from the energy gap. This leads to a convergence criterion such that once again the crossover to quadratic scaling happens when the exponent of τ Q in front of the integral is less then −2. It is easy to generalise this argument to the nth moment of the statistics of work which amounts to substituting E n p instead of E p to Eq. (A.29). As argued in the main text, this is the leading term in the nth cumulant of the distribution as well, that concludes the perturbative reasoning behind the results of Sec. 5.

B Ramp dynamics in the free fermion field theory
The non-equilibrium dynamics of the transverse field Ising chain is thoroughly studied in the literature. Due to the factorisation of the dynamics to independent fermionic modes solving the time evolution amounts to the treatment of a two-level problem parametrised by the momentum k. This two-level problem can be mapped to the famous Landau-Zener transition with momentum-dependent crossing time. Its exact solution is known and yields a particularly simple expression for the excitation probability of low-momentum modes p k (or |α(k)| 2 with the notation of adiabatic perturbation theory, cf. Sec. 2.3) in the limit τ Q → ∞. Then the KZ scaling of various quantities follows [8,13] and extends to the full counting statistics of defects [34] and excess heat [38]. For a finite Landau-Zener problem one can express the solution in terms of Weber functions [24,31] or for a generic nonlinear ramp profile as the solution of a differential equation [52,99].
To generalise the analytical solution on the chain to the free field theory we performed the scaling limit on the expressions of Ref. [52]. We remark that in the works cited above there are several parallel formulations of this problem on the chain each with a slightly different focus. Our choice to use this specific one in the continuum limit is arbitrary but the result is the same for all frameworks. We use the following notation: c k,i to refer to the operators that diagonalise the Hamiltonian initially before the ramp procedure. The operators c and η are related via the Bogoliubov transformation where the coefficients are U k = cos θ k /2 and V k = sin θ k /2 with exp(ıθ k ) = g − exp(ık) From a dynamical perspective U and V relate the adiabatic (instantaneous) free fermions and quasiparticles, hence we are going to refer to them as adiabatic coefficients. The dynamics can be solved in the Heisenberg picture using the Ansatz c k (t) = u k (t)η k,i + ıv * −k (t)η † k,i .  matrix size  25  1330  35  9615  45  56867  27  1994  37  14045  47  78951  29  3023  39  20011  49  110053  31  4476  41  28624  51  151270  33  6654  43  40353  53  207809   Table C.1: Matrix size vs. cutoff cut-off p max /m = 2π. At volume L = 50 this amounts to 100 modes in the two sectors together.
For the intensive quantities considered in Secs. 4 and 5 we worked in the thermodynamic limit L → ∞ where the sum over momentum modes is converted to an integral. Calculating the excitation probabilities of several modes up to a cutoff p max /m = 30 we used interpolation to obtain a continuous n p function. This was used in the momentum integrals that yield the energy density and its higher cumulants. The need for the higher cutoff stems from the fact that n p is multiplied with higher powers of the dispersion relation for higher cumulants.

C.1 Conventions and applying truncation
The Truncated Conformal Space Approach was developed originally by Yurov and Zamolodchikov [63,64]. It constructs the matrix elements of the Hamiltonian of a perturbed CFT in finite volume L on the conformal basis. For the Ising Field Theory the critical point is described in terms of the c = 1/2 minimal CFT and adding one of its primary fields φ as a perturbation yields the dimensionless Hamiltonian: where ∆ is the mass gap opened by the perturbation, l = ∆L the dimensionless volume parameter and ∆ φ is the sum of left and right conformal weights of the primary field φ. The matrix elements of H are calculated using the eigenstates of the conformal Hamiltonian H 0 as basis vectors: where c = 1/2 is the central charge. The truncation is imposed by the constraint that only vectors with E n < E cut are kept, where E cut is the cut-off energy. It is convenient to characterise the cut-off with the L 0 +L 0 eigenvalue N instead of the energy as it is related to the conformal descendant level. for the range of cut-offs that were used in this work. We remark that the maximal conformal descendant level N max is related to the cut-off parameter as N max = (N cut − 1)/2.

C.2 Extrapolation details
To reduce the truncation effects, we employ the cut-off extrapolation scheme developed in Ref. [73]. A detailed description of this scheme is presented in Ref. [ The exponents α < β depend on the observable O, the operator that perturbs the CFT, and on those entering the operator product expansion of the above two. For the excess energy and the magnetisation one-point function as well as the overlaps it is straightforward to apply this recipe to obtain the leading and subleading exponents. In the case of higher cumulants of the excess heat there is no existing formula. However, as they can be expressed as the sum of products of energy levels and overlaps, the leading and subleading exponents coincide with those of the first cumulant, i.e. the excess heat. The exponents are summarised in Table C.2. Sampling the dynamics using different cut-off parameters we obtained the extrapolated results by fitting the expression Eq. (C.4) to our data. In certain cases the fit with two exponents proved to be numerically unstable reflected by large residual error of the estimated fit coefficients. In these cases, only the leading exponent was used. For dynamical one-point functions the extrapolation procedure was applied in each "time slice". As evident from the exponents, the E 8 model exhibits faster convergence in terms of the cut-off. However, in most of the cases the extrapolation scheme yields satisfactory results in the FF model as well, with the notable exception of the magnetisation, as discussed in the main text. Let us now present how the extrapolation works for various quantities to illustrate its preciseness and limitations. Let us start with calculations concerning dynamics on the free fermion line. Out of the two dynamical one-point functions, the order parameter is more sensitive to the TCSA cut-off. Fig. C.1. presents an example of the cut-off extrapolation for this quantity with M i L = 50 and M i τ Q = 128. The extrapolation error (denoted by a grey band around the curve) is relatively large and partly explains the lack of dynamical scaling before the impulse regime in Fig. 4.1b. We remark that in this case the two-exponent fit was unstable hence only the leading term of Eq. (C.4) was used. The dependence on the cut-off is less drastic for shorter ramps.
The energy density exhibits much faster convergence in terms of cut-off in both models. It is in fact invisible on the scale of Figs. 4.1a and 4.2a, consequently we do not present the details of their extrapolation here. To make contrast with Fig. C.1, we illustrate with Fig. C.2 that the time evolution of the magnetisation operator is captured much more accurately by TCSA in the E 8 model. The two-exponent fit is numerically stable in this case hence we use both the leading and the subleading exponent to determine the infinite cut-off result. The change between data obtained using different cut-off parameters and the extrapolation error falls within the range of the line width in almost the whole duration of the ramp.
Apart from dynamical expectation values of local observables, we also discussed higher cumulants of work in the main text. Although the use of TCSA to directly calculate such quantities is unprecedented, based on the discussion following Eq. (C.4) we expect that the same expression accounts for the cut-off dependence as in the case of local observables. This is what we find inspecting Fig. C.3. The depicted data is a small subset of all the extrapolations whose results are presented in the main text but they convey the general message that cumulants can be obtained accurately using TCSA. The relative error in the extrapolated value is typically in the order of 1 − 3% for cumulants in the free fermion model (with an increase towards higher cumulants) and around 0.1−0.7% in the E 8 model.