Hamiltonian truncation approach to quenches in the Ising field theory

In contrast to lattice systems where powerful numerical techniques such as matrix product state based methods are available to study the non-equilibrium dynamics, the non-equilibrium behaviour of continuum systems is much harder to simulate. We demonstrate here that Hamiltonian truncation methods can be efficiently applied to this problem, by studying the quantum quench dynamics of the 1+1 dimensional Ising field theory using a truncated free fermionic space approach. After benchmarking the method with integrable quenches corresponding to changing the mass in a free Majorana fermion field theory, we study the effect of an integrability breaking perturbation by the longitudinal magnetic field. In both the ferromagnetic and paramagnetic phases of the model we find persistent oscillations with frequencies set by the low-lying particle excitations not only for small, but even for moderate size quenches. In the ferromagnetic phase these particles are the various non-perturbative confined bound states of the domain wall excitations, while in the paramagnetic phase the single magnon excitation governs the dynamics, allowing us to capture the time evolution of the magnetisation using a combination of known results from perturbation theory and form factor based methods. We point out that the dominance of low lying excitations allows for the numerical or experimental determination of the mass spectra through the study of the quench dynamics.


Introduction
In this work we investigate global quantum quenches in a quantum field theory in one spatial dimension. A global quantum quench is the non-equilibrium time evolution in isolated quantum systems which is initiated by a sudden change of one or more parameters in the Hamiltonian [1,2]. This protocol is routinely engineered in current cold-atom experiments [3][4][5][6][7][8][9][10][11], and has become a theoretical and experimental paradigm in non-equilibrium dynamics of quantum systems.
Quantum dynamics in one-dimensional integrable systems is especially interesting due to the lack of thermalisation as a consequence of the infinite number of conserved quantities present in these systems. This has been observed experimentally [3,6] and led to the theoretical proposal of the Generalised Gibbs Ensemble (GGE) [12] which has also been observed experimentally [10]. In a surprising development, recent studies [13][14][15] led to the conclusion that the completeness of the GGE requires the inclusion of novel quasi-local charges [16] which are sufficient to uniquely determine all post-quench quasi-particle distributions; indeed this latter condition had previously been shown to be necessary for the GGE to work at all [17,18].
In many cases, non-integrable systems can be considered as integrability breaking perturbations of an otherwise integrable system. Generally, non-integrable systems are expected to relax to a thermal (Gibbs) state, at least for a suitable class of (local, or few-body) observables; the principle underlying thermalisation in closed quantum systems is the Eigenstate Thermalisation Hypothesis (ETH) [19,20]. This leads to some important questions, chief among them is what aspects of integrability are still retained after such a perturbation.
A natural scenario is dubbed prethermalisation [21][22][23] which posits that for a weak breaking of integrability the system first relaxes to a deformed version of the GGE and is trapped there for some (possibly long) time before eventual thermalisation. Prethermalisation, in the sense of the system being trapped in a non-thermal metastable state, have also been observed experimentally [6]. Vestiges of integrability may also be reflected in a quantum version of the KAM (Kolmogorov-Arnold-Moser) theorem [24].
However, recent work [25] on the quantum Ising spin chain has drawn attention to the possibility that in some systems the integrability breaking interaction is never perturbative for any nonzero value of its coupling, leading to a dramatic change in the dynamics due to non-perturbative effects such as confinement.
The system considered in this work is the massive scaling Ising field theory, the scaling limit of the quantum Ising spin chain. Due to the scaling limit, quenches in the field theory correspond to very small (infinitesimal) quenches on the spin chain close to its quantum critical point. While these quenches form a somewhat special class of processes from the lattice point of view, the dynamics are still very nontrivial due to the proximity of the critical point. An important motivation for considering these quenches is that in this limit one can study aspects of quenches that only depend on the universality class of the system. Apart from this, field theoretic quenches are also interesting in their own right, since quantum field theories are valid descriptions of physical phenomena, in particular in high energy physics [26].
In the context of quantum field theory in one spatial dimension, quenches to integrable models have been widely considered [27][28][29][30][31][32][33][34][35][36][37][38][39][40]. However, much less is known about the dynamics of quenches to non-integrable quantum field theories; for integrable pre-quench dynamics, a perturbative approach was proposed in [34]. While perturbation theory has limited applicability in quantum quenches, there are some situations for which it leads to interesting predictions that can be tested within our framework.
As a result of the above situation, in this work we are particularly interested in the effects of integrability breaking in quantum field theory quenches. To address this problem, we adopt a nonperturbative Hamiltonian truncation approach that has been successfully applied to study equilibrium properties of both integrable and non-integrable two-dimensional quantum field theories.
The particular version of Hamiltonian truncation used in this work is based on a truncated fermionic space approach developed for the Ising field theory [41,42] abbreviated as TFSA. Note, however, that similar Hamiltonian truncation methods apply to a much wider range of field theory models such as perturbed minimal conformal field theories [43,44], sine-Gordon [45], Φ 4 and Landau-Ginzburg [46][47][48][49], and Wess-Zumino models [50][51][52], and can even be extended to more than one spatial dimension [53]. As a result, there are many potential directions to extend the studies in the present work.
There have also been some other recent approaches to quantum quenches incorporating Hamiltonian truncation, either combined with algebraic Bethe Ansatz methods to study one-dimensional Bose gases [24], or the truncated Hamiltonian approach forming part of a chain array matrix product state algorithm [54] to study the non-equlibrium time evolution of two-dimensional spin systems. However, our approach is essentially different from the above examples in that the simulation of the time evolution is performed within the truncated Hamiltonian setting. Since the method works with microscopic data, many quantities relevant to quench problems are accessible, such as the overlaps between states, the full statistics of work performed in the quench, the time dependent return probability (Loschmidt echo) as well as the time evolution of expectation values.
To establish the validity of this approach, first we consider the integrable case where it is possible to compare the numerical results with theoretical predictions, and then extend our studies to the non-integrable case. To have an independent verification in the latter case, the TFSA numerics are compared to infinite-volume time evolved block decimation (iTEBD) simulations on the lattice near the critical point.
The outline of the paper is the following. Section 2 contains a review of the necessary facts about the Ising field theory, and discusses the setup and general properties of quantum quenches. The Hamiltonian truncation method is introduced and described in Section 3. In Section 4, the results of the numerical time evolution applied to integrable quenches are compared to theoretical predictions. We demonstrate that it gives a correct description for integrable quenches, at least for quenches of moderate size; in particular, we avoid quenching across the phase transition. In Section 5 considers non-integrable quenches in the ferromagnetic phase. Non-integrable quenches in the paramagnetic phase are discussed in Section 6. In Section 7 we present a direct, parameter-free comparison with a lattice simulation performed for an infinite spin chain near the scaling limit. Section 8 contains our conclusions. Some of the more technical details regarding the cut-off extrapolation and finite size effects, and the description of the meson spectrum in the ferromagnetic phase are relegated to appendices.

Ising field theory and quench protocols
In this section we introduce the Ising field theory and discuss its relation to the quantum Ising spin chain. We also define the integrable and non-integrable quantum quenches that are studied in the rest of the paper.

From the quantum Ising spin chain to the scaling Ising field theory
The dynamics of the quantum Ising chain (QIC) is governed by the Hamiltonian where σ α i are Pauli matrices acting at site i, J > 0 and periodic boundary conditions σ α L+1 ≡ σ α 1 are imposed. The parameters h z and h x are called transverse and longitudinal magnetic fields, respectively.
For h x = 0 the model (2.1) simplifies to the transverse field quantum Ising spin chain that can be mapped exactly to free spinless Majorana fermions with the dispersion relation [55,56] having a gap ∆ = 2J|1 − h z |. At h z = 1 the system has a quantum critical point separating the ordered or ferromagnetic phase (h z < 1) and the disordered or paramagnetic phase (h z > 1). The order parameter is the longitudinal magnetisation σ x i that has a vanishing expectation value in the paramagnetic phase, while in the ferromagnetic phase In the scaling limit, J → ∞ and h z → 1 with fixed, the model scales to a free Majorana fermion field theory ψ(x, t),ψ(y, t) = 2πδ(x − y) . (2.5) Here we used units in which the lattice spacing is a = 2/J and consequently the resulting speed of light is c = 1. The dispersion relation of the excitations becomes (upon the substitution k = pa) that of free relativistic particles (k) → M 2 + p 2 . (2.6) We note that the distinction between the two phases is lost at the level of the field theory Hamiltonian (2.5). In fact, the two phases correspond to realisations of the same Hamiltonian on different Hilbert spaces, as described in Appendix A.1. The longitudinal magnetic field h z in the field theory couples to the spin field σ(x), which is the continuum limit of the longitudinal magnetization operator σ x i . Choosing the conformal field theory normalisation lim fixes the relation between the lattice and the continuum field as If the longitudinal field h x is non-zero, the lattice model (2.1) is non-integrable. The spectrum and the dynamics of the Hamiltonian depend strongly on the phase of the system set by h z . In the ferromagnetic phase h z < 1, the fermions correspond to domain walls, and a longitudinal magnetic fields leads to their confinement [57]. This is a non-perturbative effect which is drastic even for small values of h x , as discussed in more detail in Section 5. In the paramagnetic phase h z > 1, the elementary fermions correspond to spin waves (magnons), and a small h x only introduces some perturbative corrections to their masses and other properties.
The full scaling Ising field theory accounting for the nonzero longitudinal field has the Hamiltonian This Hamiltonian describes the universal behaviour of the spin chain (2.1) in the vicinity of the lattice quantum critical point h z = 1, h x = 0 on scales much larger than the lattice spacing a. The scaling relation between h and h x is given by wheres is defined in Eq. (2.8). For later convenience we also introduce the dimensionless combination In the ferromagnetic phase it is important to keep track of the relative orientation of the magnetic field h with respect to the spontaneous magnetisation Since in our conventions the spontaneous magnetisation is positive: 0|σ(x)|0 > 0, a field with h > 0 is in the "wrong" direction in the sense that it counteracts the spontaneous magnetisation by increasing its energy cost.

Quantum quenches in the spin chain and the field theory
In this work we consider the setting of an instantaneous global quench [1,2] which consists of an abrupt change of a parameter in the Hamiltonian. Here we study global quantum quenches in which case both the pre-quench (t < 0) Hamiltonian and post-quench (t > 0) Hamiltonian are translation invariant. As usual, the initial state at t = 0 is taken to be the ground state of the pre-quench Hamiltonian which is then evolved by the post-quench Hamiltonian, although more general choices for the initial state such as excited [58] or thermal states [59] are also possible.

Integrable and non-integrable quenches
For the quantum Ising chain (2.1), a special class of quenches corresponds to abruptly changing h z while keeping h x = 0. The dynamics of these integrable quenches can be studied analytically using free fermion techniques both on the chain [60][61][62] and in the field theory [30]. The time evolution of several quantities have been computed; we quote these results when needed in the sequel.
Another scenario that we consider is breaking integrability during the quench, in our case starting with h x = 0 and switching on a nonzero h x while also changing h z ; these quenches are called nonintegrable. In the ferromagnetic phase, integrability breaking leads to an interesting dynamics due to the onset of confinement [25]; some of the quenches considered here correspond to this case. However, the dynamics of confinement is most explicitly manifested in the evolution of the twopoint function and entanglement entropy, which is not considered in this work due to the limitations of the truncated Hamiltonian method. On the other hand, an aspect of dynamical confinement is the domination of the power spectrum by the mesonic excitations which can be studied in detail with our method thereby complementing the results of [25].
In the context of the Ising field theory Hamiltonian (2.10), the quenches we consider correspond to the parameter change (cf. Fig. 2.1) (2.13) These quenches are integrable when h = 0 and break integrability otherwise. The initial state |Ψ 0 is taken to be the ground state of the pre-quench system governed by the Hamiltonian H 0 which is H IFT in (2.10) with parameters M = M 0 and h = 0, and the time evolution is governed by the post-quench Hamiltonian H identical to H IFT in (2.10), thus (2.14) The above time evolution is numerically simulated using a truncated Hamiltonian scheme which is described in Sec. 3. Before turning to the numerical simulation, let us first recall the general expectations regarding the stationary state.

Stationary state and diagonal ensemble
After a sudden quantum quench, the system is expected to evolve towards a stationary state which is described by the so-called diagonal ensemble. Writing the initial state in terms of eigenstates |α of the post-quench Hamiltonian as the time evolution of the density matrix is given by The amplitudes C α are called overlaps and together with the post-quench energies determine the post-quench time-evolution of the system. The post-quench time evolution of an observable O is given by O(t) = Tr (ρ(t)O), while the stationary value coincides with the infinite time average (2.17) The sum (2.17) can be viewed as an ensemble average over the diagonal ensemble density matrix It is expected that a general quench results in a stationary state which is equivalent to a Gibbs ensemble for local observables O, so This is guaranteed if the so-called eigenstate thermalisation hypothesis [19,20,63] is valid. For integrable models, the stationary point cannot be described by a Gibbs ensemble [12]; the upshot of recent investigations [13-15, 64, 65] is that a suitably defined generalisation must be introduced. For integrable quenches in the Ising model, exact results for the time evolution and stationary properties have been obtained by direct calculation [30,[60][61][62]; these are used to check the validity of the TFSA in Section 4.

The truncated fermionic space approach
The main idea underlying all the variants of the Hamiltonian truncation is to split the Hamiltonian into a solvable part H 0 (usually a free theory or a conformal field theory) and a perturbation H 1 , which is assumed to be relevant in the renormalisation group sense. In order to have a discrete spectrum, the system is put in a finite volume. Using the basis spanned by eigenstates of H 0 , the matrix elements of H 1 and consequently the full Hamiltonian matrix can be computed exactly. As the spectrum of a field theory is unbounded from above, the basis must be truncated at some cut-off energy value, yielding a finite dimensional matrix. Numerical diagonalisation of this matrix gives an approximate spectrum, and expectation values or matrix elements can be computed in a straightforward manner.
The effectiveness of this method is explained by renormalisation group considerations. For a relevant perturbation the effective coupling goes to zero at high energies, and the high energy states affect only weakly the low energy properties. In any practical simulation, this can, and indeed should, be checked numerically by gradually raising the energy cut-off. Moreover, the cut-off dependence is well-understood by now through a renormalisation group in terms of the energy cut-off [53,66,67]. We emphasise that truncated Hamiltonian approaches do not rely on the integrability of the full Hamiltonian, and can be applied equally well to gapped (massive) and gapless (massless) systems.

Truncating the Hilbert space
In setting up the truncated Hamiltonian approach we follow [42] and use the eigenstates of the free fermion theory (2.5) in finite volume as a basis, truncated at a certain value of energy; this is called the truncated fermionic space approach (TFSA). Since the initial state has zero momentum and the Hamiltonian (2.10) conserves momentum, it is sufficient to include states with zero total momentum only.
Due to the finite size of the system it is necessary to specify boundary conditions. The periodic boundary conditions on the chain allow for periodic or anti-periodic boundary conditions for the Majorana fermion field, splitting the Hilbert space into the Ramond (periodic) and the Neveu-Schwarz (anti-periodic) sectors (see Appendix A.1).
The eigenstates of the free fermion Hamiltonian (2.5) can be constructed using the mode expansion where a(θ n ), a † (θ n ) = δ n,n . (3. 2) The momentum and energy of the quasi-particles are quantised as where n is the momentum quantum number and we also introduced the rapidity parameter θ for later convenience. The quantum number n is integer in the Ramond sector and half-integer in the Neveu-Schwarz sector, and we use the phase factors ω = e iπ/4 andω = e −iπ/4 following the conventions of [42]. All physical quantities are measured in units of the post-quench fermion mass M : energies are simply given in units of M while volume is parameterised by the dimensionless variable M L.
Using the mode expansion (3.1) and the rapidity variables in (3.3), the free fermion Hamiltonian (2.5) can be rewritten as The eigenvectors of this Hamiltonian are given as with energy eigenvalues M n cosh(θ n ), and the vacuum defined by a(θ)|0 = 0. The truncation of the Hilbert space is performed by keeping states with energies less than an imposed cut-off Λ: The number of states retained is denoted by N cut . For the construction of the matrix elements of the Hamiltonian (2.10) and the calculation of expectation values we need the finite volume matrix elements of the energy operator ε(x) = i :ψ(x)ψ(x) : (3.7) and the magnetisation operator σ(x). The matrix elements of ε(x) can be computed in a straightforward way from the mode expansion, while those of σ(x) were given in [42]; for the sake of completeness and in order to fix our conventions they are summarised in Appendix A.2.

Choice of the TFSA basis
In order to improve the accuracy of the TFSA, one can optimise the choice of basis when implementing the quench from H(M 0 , h = 0) to H(M, h) so that the initial state |Ψ 0 and the post-quench Hamiltonian H are both represented with the highest possible precision. This is best achieved by using the eigenbasis of the Hamiltonian H(M, 0), i.e. the free fermion theory closest to the actual post-quench Hamiltonian (see Fig. 2.1). For integrable quenches from H(M 0 , 0) to H(M, 0) the above choice of basis implies that the time evolution operator e −iH(M,0)t is diagonal and exact apart from the truncation. The overlaps between the ground state of H(M 0 , 0) with the eigenstates of H(M, 0) are exactly known (cf. Eq. (4.6)). The accuracy of the TFSA can be verified by numerically computing the ground state of the pre-quench Hamiltonian written as in the post-quench basis and comparing its overlaps with the basis states to the exact predictions.
As discussed below, we found that the TFSA approach reproduces these overlaps to a good accuracy for quenches of moderate size.

TFSA implementation of the time evolution
For integrable quenches to H(M, 0), the time evolution operator e −iH(M,0)t is diagonal in the eigenbasis of H(M, 0) we use and can be trivially calculated. This is not the case for non-integrable quenches to H(M, h = 0). For these quenches the action of e −iH(M,h)t on the initial state can be evaluated efficiently using an expansion in terms of Chebyshev polynomials, which is known to give the best finite order polynomial approximation of continuous functions on a finite interval. More details are given in Appendix A.3. As we are are principally interested in the infinite volume behaviour, the simulation time is bounded by t = L/2, i.e. by the time needed for the fastest excitations moving along the light cone to travel half-way around the circle and meet again. For t > L/2 the behaviour of observables is affected by partial revivals absent in the infinite volume case. Note that numerical errors accumulated during the TFSA simulation may in principle invalidate the results before this time is reached.

Integrable quenches
Before applying the TFSA to non-integrable quenches, we first study integrable quenches from H(M 0 , 0) to H(M, 0). As discussed in Section 2.2.1, there is a large number of analytic results for these quenches, therefore they can serve as a benchmark of the TFSA simulation. Following the discussion in Section 3.2, we work in the basis of eigenstates of the post-quench free fermion Hamiltonian H(M, 0) and find the initial state as the ground state of the pre-quench Hamiltonian H(M 0 , 0) that is numerically computed in this basis. The time evolution is trivial to implement since e −iH(M,0)t is diagonal and its eigenvalues are exactly known. We test the TFSA by comparing    its results to analytic expressions for the overlaps and for the post-quench dynamics of the energy operator, the magnetisation, and the Loschmidt echo.

Time evolution of the energy operator
We start our investigation of the out of equilibrium dynamics of various observables by studying the so-called energy operator, ε(x) = i :ψ(x)ψ(x) :, which is the continuum limit of the transverse magnetisation on the lattice. The known analytical result on the lattice [62] has the scaling limit ε(t) = Ψ(t)|ε|Ψ(t) = − n (cos β n cos ∆ n + sin β n sin ∆ n cos(2E n t)) , with η = ±1 in the paramagnetic and ferromagnetic phase, respectively. In (4.1), the summation over momentum quantum numbers is logarithmically divergent but the expression can be evaluated by imposing a momentum cut-off. In order to compare the TFSA simulations with this analytical formula, we regularise the expression by introducing a cut-off n cut for the momentum quantum numbers such that Here the factor of 2 reflects the fact that the initial state consists of particle pairs. This is a sort of a "single particle" cut-off and as such it does not exactly correspond to the "many-body" TFSA truncation specified in (3.6). In spite of this mismatch, for quenches of moderate size a very good   agreement is observed, as shown in Fig. 4.1a. For large quenches, as shown in Fig. 4.1b, the analytical formula yields high frequency oscillations that are present in the result of the TFSA simulation, but with a much smaller amplitude. However, apart from this mismatch the agreement with the exact result is still quite satisfactory.
It is immediately apparent from Fig. 4.1 that the frequency of the high-frequency oscillations is related to the cut-off. Numerical checks reveal that their frequencies are indeed very close to Λ. These oscillations are present in the numerical TFSA results obtained for ε(t) and the Loschmidt echo in the integrable quenches; for the non-integrable quenches some remnants of them survive depending on the dynamics as discussed later. The cut-off related oscillations can be suppressed by the Gaussian smoothening procedure described in Appendix B. For the energy operator, however, we do not apply smoothening as the theoretical prediction itself involves a cut-off due to the logarithmic ultraviolet divergence and consequently it also features these oscillations. The logarithmic divergence is also manifest in the data shown in Fig. 4.1: apart from the cut-off related oscillations, ε(t) attains a stationary value for M t 5, but this value steadily increases with the cut-off. Due to this divergence it makes no sense to perform a cut-off extrapolation in this case, unlike for other observables.

Statistics of work and the Lochsmidt echo
The statistics of work is an important characteristic of quantum quenches which gives the distribution of work performed during the quench [68]. For a system with discrete energy levels, the probability of performing work E α − E 0 is where the initial state |Ψ 0 is the ground state of the pre-quench Hamiltonian, and |α are the eigenstates of the post-quench Hamiltonian with energies E n . In other words, the statistics of work is given as the set of pairs (E n , P n ) of possible energy values and their respective probabilites P n = | Ψ 0 |n | 2 .
The overlaps appearing in (4.4) are exactly known for integrable mass quenches (M 0 , 0) → (M, 0), which are continuum limits of quenches of the transverse magnetic field in the Ising quantum spin chain. The components of the initial state in the Ramond and Neveu-Schwarz sectors can be expanded in terms of the Fock eigenvectors (3.5) of the post-quench Hamiltonian as [30,61,62] where the Ramond component only exists for quenches starting in the ferromagnetic phase and the normalisation factors are given by The initial state is thus a combination of zero-momentum particle pairs. In the paramagnetic phase the ground state only has components in the Neveu-Schwarz sector. In the ferromagnetic phase in infinite volume the ground state is doubly degenerate, but in finite volume this degeneracy is split and the eigenstates of the Hamiltonian have zero magnetisation. In order to study the time evolution starting from a state with finite magnetisation, the initial state must be chosen as a linear combination of the lowest energy states in the Ramond and Neveu-Schwarz sectors with equal weights: The reliability of the TFSA can be checked by comparing the numerically obtained work statistics to the analytical result following from expansion (4.6) in Fig. 4.2. Within the volume and cut-off range considered, the ground state always has the largest overlap with the initial state, while the arcs in the figure correspond to the 2-particle, 4-particle etc. states. The TFSA is found to reproduce the statistics of work well for moderate quenches. Moreover, raising the energy cut-off improves the accuracy of the overlaps given by the TFSA. Note that even though the ratios of the overlaps were perfectly reproduced by the TFSA, the truncation of the Hilbert space changes the overall normalisation of the overlaps compared to the exact prediction due to the omitted states.
A further check is provided by the Loschmidt echo or fidelity defined as Physically, L(t) is the probability of observing the system in the initial state at time t after the quench. It is also the modulus square of the characteristic function of the work statistics (4.4) (4.10)  Loschmidt echo After the cut-off extapolation, the TFSA result for the time-dependent Loschmidt echo agrees well with the theoretical prediction for moderate mass quenches within the same phase. In addition, the thin dashed line in Fig. 4.3b shows the TFSA result obtained by using the exact overlaps computed from Eq. (4.6). This eliminates completely the numerical error from the overlaps, but still leaves the error coming from truncating the basis to a finite number of vectors. For the case of non-integrable quenches considered later this leads to the important observation that using the exact overlaps (as explained in Sec. 3.2) improves the accuracy considerably, since it takes into account the integrable part of the quench almost exactly.   the order parameter initially increases, while for quenches towards the critical point it decreases. In both cases, after an initial oscillating transient an exponential decay sets in. This is in agreement with the analytical results both on the lattice [61] and in the Ising field theory [30]. In the latter, for small quenches a resummed form factor expansion in Ref. [30] yielded

Time evolution of the order parameter
where K(θ) is the kernel function in the expansion of the initial state (4.6) and α was conjectured . This expression agrees with the continuum limit of the lattice result, and apart from the subleading prefactor λ(t) it can also be obtained using semiclassical arguments [69][70][71][72]. The analytic result (4.12) (without λ(t)) is plotted using green dots in Fig. 4.4. The agreement with the extrapolated TFSA data is quite good. The decay rates τ −1 TFSA and the amplitudes A TFSA can be extracted fitting our results for 8 < M t < 20 by an exponential A TFSA exp(−t/τ TFSA ), and they are in good agreement with the exact values τ −1 andσ for both quenches: Similarly to the Loschmidt echo, Fig. 4.4b also shows in thin dashed line the TFSA result obtained using the exact overlaps that agrees almost perfectly with the analytic result.  A time dependent oscillatory prefactor can also be extracted from the TFSA results by dividing it by the asymptotic expressionσe −t/τ . As can be seen in Fig. 4.5, this prefactor agrees in its main features with with the prediction λ(t) in Eq. (4.14), in particular, its Fourier analysis reveals a single strong peak at frequency 2M. The high frequency oscillation visible in the TFSA results is a cut-off effect.

Summary of results for integrable quenches
To sum up, we have demonstrated that for quenches of moderate size within the same phase of the Ising field theory, the TFSA method is able to reproduce the theoretical results for various quantities, including the statistics of work P (W ), the Loschmidt echo L(t), and the expectation values ε(t) , σ(t) to a good accuracy.
We emphasise that good agreement was reached between the theoretical predictions valid in infinite volume and our finite volume TFSA calculations performed at system size M L = 40. In particular, this implies that the cut-off extrapolated numerical results have small finite size effects, which is expected for massive field theories that typically show exponentially suppressed finite size effects, as noted previously in [54].
Finally we note that the TFSA is less accurate for larger changes in the Hamiltonian parameters, in particular, reliable simulations for quenches across the critical point featuring non-analyticities in the Loschmidt echo would require much larger energy cut-offs.

Non-integrable quenches in the ferromagnetic phase
Having tested the TFSA method let us turn now to the case of non-integrable quenches. As discussed in Sec. 3.2, for quenches from H(M 0 , h 0 = 0) to H(M, h) it is possible to use the eigenbasis of the Hamiltonian H(M, 0) and construct the initial state truncating the expressions (4. 6,4.8).
In this Section we consider quenches within the ferromagnetic phase where the time evolution is found to be dominated by the meson bound states. They appear in the spectrum as a result of turning on the external magnetic field in the longitudinal direction.  In the ferromagnetic case seven mesons are visible below the continuum, the blue dashed lines show the first three meson masses computed for an infinite system. The state with linearly increasing energy is the false vacuum corresponding to the ferromagnetic ground state with negative magnetisation. In the paramagnetic case there is a singe massive particle.

Spectrum
Before turning to the quenches, we discuss the spectrum of the post-quench Hamiltonian. In the ferromagnetic phase, switching on a longitudinal magnetic fieldh causes a non-perturbative change in the spectrum which can be understood most easily on the lattice. In the ferromagnetic phase, the elementary excitations ath = 0 corresponding to free fermions are domain walls. A longitudinal field induces a linear potential between pairs of domain walls that are separated by a distance d, which for a small longitudinal magnetic field h has the form leading to a confinement of the fermions into so-called mesons, a scenario first proposed by McCoy and Wu [57]. Based on this picture, it is possible to compute the spectrum using various theoretical approaches, for details see Appendix D.
The spectrum of the first few zero-momentum excited states as a function of the dimensionless system size M L is shown in Fig. 4.6a. Note that the energy values quickly reach a stable value, already around M L = 30, which are close to the infinity volume values even at this relatively small cut-off. The plot shows the values of the first three meson masses for comparison as computed using the WKB approximation (see Appendix D). The excellent agreement demonstrates both the validity of the WKB approximation and the precision of the TFSA method.

Time evolution of the order parameter
Turning to the time evolution of the expectation value of the magnetisation operator σ after the quench, we present numerical data for quenches from M 0 = 1.5M to M while the magnetic field is changed fromh 0 = 0 to eitherh = −0.05 (Fig. 5.1a) orh = −0.1 (Fig. 5.1b). A larger value ofh causes the features of integrability breaking to appear on a shorter time scale. 1 The extrapolation introduces spurious oscillations at short times that are not present in the raw data. Smoothing the raw curves using a Gaussian kernel before extrapolation (see Appendix B.2) eliminates these oscillations (see blue dots in Figs. 5.1a,b) but spoils the large time behaviour.
Already a small magnetic field results in a drastic change in the time evolution of the magnetisation. Instead of the exponential decay seen for integrable quenches withh = 0, it does not decay to zero but features slow and fast oscillations. The various oscillation frequencies can be extracted by taking the Fourier transform of the time series σ(t) in a time window 4 < M t < 36 in order to avoid the spurious oscillations at short times induced by the extrapolation procedure. As demonstrated in Figs. 5.1c, 5.1d the resulting power spectrum shows clear peaks related to the meson masses, suggesting that the confinement has drastic effects on the quench dynamics. Note that a non-zero asymptotic magnetisation is expected due to the presence of the external magnetic field. Both the overall behaviour of the magnetisation and the relation between the oscillation frequencies and the meson masses are in full agreement with the results found in the earlier work [25] which considered the quench dynamics on the spin chain governed by the Hamiltonian (2.1).
Fig. 5.2 shows σ(t) computed in a larger volume M L = 100 so that the time evolution can be followed further, up to M t = M L/2 = 50. The magnetisation is seen to oscillate around its infinite time average given by the diagonal ensemble average. The oscillations persist up to the latest times we were able to study, even when the system size L is much bigger than the equilibrium correlation length M −1 . Note that the studied times are already long in the natural unit M −1 , so any qualitative change will supposedly take place at some much larger, possibly astronomical time scales, if ever. The questions whether these oscillations eventually die out and if they do, what is the time scale of this damping, should be the subject of further study.  The blue dots at short times are obtained by smoothening the data before extrapolation.

Statistics of work and Loschmidt echo
The dominant role played by mesons is further supported by the statistics of work which is shown in Fig. 5.3a for a specific but illustrative quench in the ferromagnetic phase. It is clear that a few of the lowest lying states in the spectrum dominate and are responsible for the oscillatory behaviour, while all other states have much smaller overlap with the initial state.
The dominance of low lying states runs contrary to normal expectations: in a quantum quench usually the highly excited states in the middle of the post-quench spectrum are the relevant ones. On the one hand, the agreement with infinite volume lattice simulations indicates that this feature is robust against changing the cut-off or the volume of the system. On the other hand, as discussed in detail in Appendix C, the overlaps of the 1-meson states |Ψ m L are eventually expected to decay as for large enough L. Although in our simulations we cannot access such large volumes, the theoretical curve fits very well the overlaps measured by TFSA (see Fig. 5.4). However, as can be seen from the table in Fig. 5.4, the critical volume L crit = 1/B at which the one-particle overlap is expected to start decreasing is a few hundred times the correlation length. In such large volumes the finite size effects are already expected to be minuscule; eventually, as shown by the excellent agreement with the iTEBD lattice simulations discussed in Section 7, already the volume M L = 40 which is just about a hundred times the correlation length can be considered infinite in terms of field theoretic finite size effects. This can be explained by observing that the smallness of the exponent B, which from Eq. (C.14) is given by which is small whenever the post-quench density of particles given by [61] ρ =ˆ∞ is small. For a small density of particles, the average distance between particles is much larger than the correlation length, so even though the infinite volume post-quench contains a finite density of particles, its dynamics is dominated by one-particle contributions "semiclassically pasted together" [71][72][73]. This opens the way to measure the masses of the mesons through studying the post-quench time evolution of the magnetisation. This method was called "quench spectroscopy" in Ref. [74]. Fig. 5.3b shows the Loschmidt echo obtained by the TFSA. It quickly relaxes from 1 to some lower value and then begins to oscillate with an amplitude that increases in time. Note that the extrapolation induces fast oscillations for short times that are most likely numerical artefacts. These can be removed by Gaussian smoothening, as plotted in blue dots in Fig. 5.3b, which however spoils the steep decreasing edge of the initial transient. In any case, the curves seem to be reliable after t ≈ 3.

Non-integrable quenches in the paramagnetic phase
In this section we review the results obtained for quenches within the paramagnetic phase. This phase does not feature confinement and consequently has a very different excitation spectrum containing magnons and no mesons. As a result, the time evolution of the magnetisation and the Loschmidt echo significantly differs from those in the ferromagnetic phase.

Spectrum
Similarly to the ferromagnetic case, we first show the spectrum of the post-quench Hamiltonian in Fig. 4.6b. The first excited state is the magnonic wave with a mass perturbatively corrected by the magnetic field. It approaches its infinite value as fast as the meson states in the ferromagnetic phase (cf. Fig. 4.6a). The higher excited states shown in the plot correspond to two-particle states of the magnons and form a continuum in the large volume limit. There are also states with more than two magnons that are not shown in the plot.

Time evolution of the order parameter
To make contrast with the time evolution of σ(t) in the ferromagnetic phase that was dominated by the presence of a few low energy meson bound states arising due to the confinement of kink-like excitations, we examine the dynamics of σ(t) in the paramagnetic case where no such mesons are present. As shown in Fig. 6.1, σ(t) shows a slightly damped sinusoidal oscillation.
In contrast to the ferromagnetic phase where the dynamics are dominated by the non-perturbative effects of confinement, in the paramagnetic phase a perturbative treatment of the post-quench dynamics is expected to give a reasonable description of small quenches in h (without changing the mass, i.e. M 0 = M ). Such an approach was worked out in Ref. [34]. We take the pre-quench Hamiltonian H(M, 0) to be the unperturbed Hamiltonian and ∆H = h´dx σ as the perturbation. Using the results of [34], the leading perturbative contribution to the time evolution of the magnetisation is given by the one-particle term where F 1,0 = p = 0|σ|0 is the infinite volume form factor of σ between the ground state |0 and a zero-momentum one-particle state taken in the pre-quench basis (c.f. (A.9)). This formula does not account for the constant offset that is however expected due to the non-zero magnetic field.  Table 6.1: Comparing fitted parameters of the σ(t) curves in the paramagnetic phase with other TFSA and perturbative data. σ DE is the diagonal ensemble prediction for a mass quench M 0 = 1.5M at cut-off Λ/M = 8. A and f are fitted using the data for a quench with M 0 = M while C is fitted for a quench with M 0 = 1.5M .
Moreover, for simultaneous quenches both in h and M the oscillation is seen to be damped (c.f. Fig.  6.1). A natural function to fit the TFSA result is For quenches in h with M 0 = M there is no damping, so τ = 0. In Table 6.1 we report the fitted values of the amplitude A together with the perturbative prediction in (6.1), the frequency ω with the h-dependent mass of the magnon, and the offset C along with the expectation value of σ in the diagonal ensemble (infinite time average). The agreement for all three parameters is excellent, showing that the offset is given by the infinite time average, and perturbation theory captures both the amplitude and the frequency of the oscillations although for the latter the post-quench value of the particle mass gives a better agreement.
The mass correction can also be estimated in the framework of second order form factor perturbation theory [75] (including intermediate states with up to two particles) with the result which is quite close to the value extracted from fitting the mass gap data in Table 6.1 up toh = −0.05: from ω fit . (6.5) This shows that the quench spectroscopy works also in the paramagnetic phase. Apart from the mass of the particle excitation, the quench dynamics also allows for the determination of the form factor appearing in the amplitude in (6.1). For combined quenches both in h and M, we found that the function (6.2) with τ > 0 describes almost perfectly the evolution of the magnetisation, at least for moderate size quenches (see Fig.  6.1). The frequency of the oscillation is found to be largely independent of M 0 , it changes around 1% while M 0 is varied between 0.5M and 1.9M. The change in the offset C is around 4 − 5%, while the amplitude is a bit more sensitive to the initial mass.
The damping rate, τ −1 , is shown in Fig. 6.2 for various values of h and M 0 . We found that it depends only very weakly on the value of h for small h. Moreover, it is close to the damping rate valid for integrable mass quenches at h = 0, given by Eq. (4.13) and shown in dashed line in Fig. 6.2. For integrable quenches within the paramagnetic phase the magnetisation is zero for all times, but expression (4.13) holds for quenches from the ferromagnetic phase to the paramagnetic phase [61]. Interestingly, we found that the same expression gives the decay rate for quenches within the paramagnetic phase when the post-quench magnetic field is non-zero. The reason for this is that the main source of damping is the finite density of magnons which is predominantly induced by the quench in the mass, described by Eq. (4.13). The small deviations are due to higher order effects both in h and in ∆M.

Statistics of work and Loschmidt echo
As we saw above, the time evolution of σ is dominated by a single frequency in the paramagnetic case. To support this observation, in Fig. 6.3a we show the statistics of work in this phase that is similar to the ferromagnetic case in that it is dominated by low lying states. The difference is that there is a single massive particle leading to a pure oscillation instead of multiple frequencies. In Fig.  6.3b the Loschmidt echo is plotted in the paramagnetic phase featuring a persistent oscillation with a single frequency.

Comparison with iTEBD lattice simulations
As a final check of our TFSA method, in this section we compare our results for the time evolution of the magnetisation with simulations on the Ising spin chain with both transverse and longitudinal fields governed by the Hamiltonian (2.1). On the lattice we used the infinite volume Time Evolving Block Decimation (iTEBD) method which is free of finite size effects and thus allows for testing the finite size errors of our TFSA results.
In order to carry out the comparison, the lattice simulation should approach the scaling limit, i.e. it must be performed at sufficiently large values of J near the quantum critical point such that  In Fig. 7.1a the results for σ x n (t) of three iTEBD simulations are shown with J = 5, 10, 20 rescaled according to Eq. (7.1) such that they all correspond to the quench H(M 0 = 1.5,h = 0) → H(M = 1,h = 0.1) in the ferromagnetic phase. The three curves are almost indistinguishable from each other, so they are very close to the universal curve corresponding to the field theory. In Fig.  7.1b we compare this curve with our TFSA result and find excellent agreement. We stress that this comparison has no free parameters. Similar comparisons are shown for a ferromagnetic quench in the h direction without a mass quench in Fig. 7.1c and for a combined quench both in M and h in the paramagnetic phase in Fig. 7.1d. For the quench only in h shown in Fig. 7.1c some deviations can be seen. However, the magnitude of these deviations is very small (note the scale on the vertical axis). Moreover, the discrepancy is largely eliminated by modifying relation (7.1) by a factor ≈ 0.9998, the corresponding curve is shown in thin black line.
From the very good agreement between the lattice and the TFSA results we can draw the following conclusions. First, as the iTEBD does not suffer from finite size corrections, the finite size errors in the TFSA are very small, presumably exponentially suppressed. The finite size L however restricts the simulation to times t ≤ L/2. Secondly, since the iTEBD method truncates the Hilbert space based on entanglement entropy instead of energy and this truncation is very well controlled, we conclude that the extrapolation in the TFSA cut-off has a rather small error and eliminates the effect of the finite energy cut-off to a great extent.     Let us emphasise that the TFSA is much more efficient to study the dynamics of the field theory than the iTEBD on the lattice in the scaling limit. The reason is that approaching the quantum critical point, the lattice simulation takes more and more time due to the critical slowing down. In practice, this means that each simulation of a given quench can take days or even weeks depending on the distance from the critical point (i.e. the value of J) on an Intel 3.60 GHz processor. In contrast, in the TFSA the most time consuming step is the computation of the matrices that can take days for Hilbert spaces of dimension ∼ 10 5 . However, once the matrices are constructed they can be used to investigate many different quenches as the simulation of the time evolution only takes a few hours.

Conclusions
In this paper we studied the real time evolution following quantum quenches in the scaling Ising field theory using the Truncated Hilbert Space Approach. Our main results are twofold.
On the one hand, we demonstrated that the TFSA method is suitable for studying the quench dynamics in one dimensional interacting field theories. Due to the microscopic nature of the method it has access to many quantities, including the individual overlaps, the work statistics, or the Loschmidt echo that are hard to access with other methods. As there is no universal numerical method available for studying the out of equilibrium dynamics in continuum systems, the present proof of principle study can have far reaching consequences. This is especially so since using already existing truncated Hamiltonian approaches, the method can be easily generalised to other field theories. An obvious candidate is the sine-Gordon model for which the method was first developed in [45] and which is an important model providing the low energy effective description of various one dimensional systems, including coupled 1D quasi-condensates currently studied in matter wave interference experiments [76]. The method can also be generalised to arbitrary finite duration quench protocols or ramps.
On the other hand, we found interesting dynamics in the non-integrable Ising field theory. For quenches in the ferromagnetic phase we found that the weak confinement of the elementary excitations of the integrable model has drastic consequences for the quench dynamics, in accordance with a recent study of quenches in the lattice version of the model [25]. We also found that in both phases the quench dynamics is dominated by a few low energy states even for relatively long times. In the paramagnetic phase this can be explained borrowing a perturbative approach developed by Delfino in [34], while in the ferromagnetic phase this turns out to be connected to low density during the quench, leading to the one-particle states dominating the system up to system sizes of hundreds of the correlation length. This means that measuring specific observables either in numerical simulations or in experiments allows one to extract important information about the spectrum and form factors of the post-quench Hamiltonian, realising a kind of quench spectroscopy.
Finally, the agreement between the lattice iTEBD and the continuum TFSA also shows that the scaling limit of quantum quenches on the lattice agrees with sudden quenches in the field theory. Therefore, the universal description given by the scaling field theory in the vicinity of the critical point is valid beyond the equilibrium and captures the non-equilibrium dynamics of the spin chain.

A Some details of the finite volume Ising field theory
In this appendix we review some basic ingredients of the finite volume Ising field theory, such as the structure of the Hilbert space in the two phases, and the finite volume matrix elements of the magnetisation operator. We also provide some details on implementing the time evolution using the Chebyshev expansion.

A.1 The ferro-/paramagnetic phases in field theory
As discussed in the main text, the free fermion Hamiltonian (2.5) is the same in both the paramagnetic and ferromagnetic phases. To understand the eventual difference between the phases at the level of the field theory, one can consider the Jordan-Wigner transformation from the spin chain (2.1) to the fermionic formulation. This maps the Z 2 symmetry of the Ising chain into the fermionic parity e iπN , where N is the fermion number operator. The mapping to fermions specifies the boundary conditions to be periodic in the odd N sector and anti-periodic in the even N sector. Consequently, the set of possible wave numbers is given by integer multiples of 2π/L for N odd, and half-integer multiples for N even, corresponding to the Ramond (R) and Neveu-Schwarz (NS) sectors: The excitation energy of the lowest lying excitation compared to the Fock vacuum is proportional to h z − 1 which becomes negative in the ferromagnetic phase, so the ground state must be redefined by filling this negative energy level. As a consequence, the Ramond sector in the ferromagnetic phase has states with even number of particles, in contrast to the the paramagnetic phase. The number of particles in the Neveu-Schwarz sector is even in both phases. To summarise, the basis used in the ferromagnetic phase is given by while the basis in the paramagnetic phase is where |0 R and |0 NS are the Ramond and Neveu-Schwarz vacua, respectively, and a † p creates a fermion particle with momentum p and mass M .

A.2 Matrix elements of σ
In the Ising field theory, the order parameter field σ only has matrix elements between the even parity (Neveu-Schwarz) and odd parity (Ramond) sectors. In the eigenbasis of H(M, 0) in Eq. (2.5), the matrix elements have the following exact expression [42,77,78]: where θ n is the rapidity corresponding to momentum p = 2π/L as defined in (3.3b). F m,n is the infinite volume matrix element (form factor) given by F m,n (θ 1 , . . . , θ m |θ 1 , . . . , θ n ) with [x] being the integer part of x, andσ is the infinite volume expectation value of the order parameter field given in (2.9).

A.3 Implementing the quench time evolution
The Chebyshev polynomials T n (x) are defined by the recurrence relation (A.14) They satisfy the orthogonality condition and form a complete basis for functions on the interval [−1, 1]. This can be used to evaluate functions of matrices whose eigenvalues lie within the unit circle. In particular, the exponential function appearing in the time evolution operator can be expressed as where J n (z) are the Bessel functions (A. 17) In order for the above expansion to be valid the Hamiltonian H must be shifted and rescaled so that all its eigenvalues lie inside the unit circle; the effect of this transformation on the time evolution can be undone by an appropriate rescaling of the time t. In a practical calculation the expansion is truncated at some n; for our purposes we used a truncation at n = 1000.
In practice, there is no need to evaluate the time evolution operator itself, only its action on the initial state, that is, the time evolved wave function. Thanks to the recursion relations (A.14), this only requires matrix-vector multiplication instead of matrix-matrix multiplications.
Once we have the time evolved wave function at a given time as an N cut dimensional vector, the calculation of expectation values simply amounts to multiplication between this vector and the numerical matrix of the operator.

B Cut-off extrapolation and finite size effects
In this appendix we discuss the cut-off effects in the TFSA method and provide details on the extrapolation in the cut-off.

B.1 Renormalisation group improvement of TFSA
Since our interest is in the continuum field theory, it is desirable to eliminate the cut-off dependence inherent in the TFSA. This is achieved by introducing running parameters m and h that are adjusted according to suitable renormalisation group equations; this idea was recently introduced in various truncated Hamiltonian approaches to field theories [53,66,67].
The Ising field theory can be represented as a relevant perturbation of its ultraviolet fixed point: where the fixed point Hamiltonian is the conformal field theory of the massless Majorana fermion, and The leading order renormalisation group equations for the couplings τ and h are [79] dh dn =hτ 2π n 2 , dτ dn =h To the leading order, the TFSA has the same running for the couplings since the introduction of mass affects only low-lying states and is a sub-leading effect, therefore we only need to transcribe the above RG equations in terms of the parameters of the TFSA. The cut-off in the free fermion theory can be written in terms of the conformal cut-off as Λ = 4πn/L; in order to have a dimensionless cut-off λ and system length l we introduce some mass scale m (to be specified later): Using these relations we can write the effective couplings at some value of the free fermionic cut-off λ in terms ofτ so in order to ensure m = M (λ) the value of l must be adjusted as a function of the cut-off. This is very inconvenient as it means that the basis for the TFSA must be regenerated each time the cut-off is adjusted, which is quite time consuming. On the other hand, choosing l to be independent of λ, the RG equations can be rewritten in terms of the dimensionless parametersm  The numerical solution of the above equations shows that the parameter h changes significantly while flowing from λ = ∞ to the values of the cut-off used for the TFSA simulations (λ < 10). However, to a good approximation M remains constant which means that using M = 1 without adjusting l is a very good approximation. This leads to the choice m = M , (B.10) i.e. we choose units in which M = 1.
In contrast, the value of h must be adjusted according to the RG flow; whenever we refer to the value of h in the main text it means the physical value at infinite cut-off (λ = ∞) while the numerical simulations use the value corresponding to the cut-off implemented in the TFSA.

B.2 Cut-off dependence and extrapolation schemes
Since the cut-off is imposed in the total energy of the states in the TFSA space, the dimension of the Hilbert space increases quickly with the volume L. As a consequence the maximum value of Λ allowed by the available computing power is smaller for larger systems and the actual maximum values of Λ for different values of L are given in Table B.1. We made an exception for the overlaps with the initial state, where for L = 30 the maximum cut-off was extended to Λ = 11 (5.3a).
To eliminate the cut-off dependence remaining after the introduction of running couplings, we performed an extrapolation of certain quantities in the cut-off. For this purpose we used a series of cut-offs up to the maximum values indicated, with a spacing of ∆Λ = 0.25. The extrapolation was always performed using simulations with the 6 largest available cut-offs, for example in Fig. 5.1a the cut-offs Λ = 8.5, 8.75, . . . , 9.75 were used.
The extrapolation itself used a certain fitting function for the cut-off dependence. For the evolution order parameter σ(t) there is a theoretical prediction available from the formalism introduced in [80], which predicts that the leading cut-off dependence is given by   gives a very accurate fit to the cut-off dependence, as illustrated in Fig. B.3. However, before extrapolation the oscillations related to the cut-off must be eliminated, since their frequency depends on Λ which leads to some unphysical features after extrapolation. We eliminate them by convolving the curves σ(t) Λ , L(t) Λ with a Gaussian where the parameter Σ is adjusted according to the cut-off, Σ = 2π/Λ. Note that at the edges of the time window the Gaussian function would involve a substantial contribution from nonexistent data points that are outside the window; this is avoided by adjusting the width Σ for time instances close to the edges of the window.

C Finite size behaviour of the overlaps
Here we present numerical data for the meson wave functions in order to show that two-quark states dominate the meson states. Based on this observation we present a theoretical calculation of the finite size dependence of the overlaps of the 1-meson states. where x is the relative coordinates of the two quarks. Standard quantum mechanics then yields the normalised wave function where Ai is the Airy function, with z > 0 corresponding to a zero of the Airy function There is an infinite number of zeros which can be ordered as z 1 < z 2 < . . . , with z n giving the nth meson state. The variable E gives the mass m of the meson via Introducing the momentum space wave function (with a phase redefinition to make the result purely real) the properly normalised finite volume wave function can be written as Note that the index here runs over both integers and half-integers, corresponding to the Ramond and Neveu-Schwarz components of the state, respectively. Using the TFSA, the meson wave function can be explicitly evaluated in the fermionic basis. We have performed this calculation for values of longitudinal magnetic fieldh = −0.05, −0.1, and volumes M L = 30, 40. The first important observation is that in each case, the two-quark components accounted for more than 99% of the squared norm, so the two-quark approximation is very accurate.
The finite volume wave function can be directly compared to the numerical components of the meson state vectors computed in the TFSA basis via the relation Ψ(q) = √ LΨ L (q) (C.9) As illustrated in Fig. C.1, we found good agreement which gets worse for the higher mesons (and also with increasingh), primarily because the non-relativistic approximation is less accurate. Note also that the finite volume wave functions corresponding to different volumes scale on top of each other, and the momentum range where the wave-functions have their support lies well below the cut-off, which excludes the presence of substantial finite size or volume dependence remaining. Therefore the deviations of the TFSA wave function from the prediction Eq. (C.7) observed in which is much larger than the volumes that can be reasonably reached using TFSA, and is also too large to have any finite size dependence for field theoretic correlators. The consequences of this observation are discussed in the main text.

(D.3)
Another approach to compute the meson mass is by solving a Bethe-Salpeter equation [81][82][83], but it is only accurate for weak magnetic field h, while the semiclassical results give a very good approximation for a wide range of parameters as demonstrated in [84].

E Some details of the iTEBD simulations
In order to obtain an external confirmation of our calculations with independent methodology we employed iTEBD (infinite Time Evolving Block Decimation) numerical simulations to follow the entire relaxation process toward the equilibrium state. The algorithm is based on the infinite Matrix Product State (iMPS) description of one-dimensional translational invariant lattice models in the thermodynamic limit which is free of any finite size effect. The canonical iMPS representation of a generic many-body state is therefore Λ e · · · ] |· · · s j s j+1 · · · , (E.1) where Γ s j o/e are χ × χ matrices associated with odd/even lattice sites, with s j spanning the Hilbert space of the j th site with canonical basis {|↑ , |↓ }; similarly, Λ o/e are diagonal matrices with entries being the singular values associated with Schmidt decomposition for the bipartition of the system on the odd/even bonds. Notice that the iMPS representation is properly defined thanks to the right/left orthonormalisation conditions where h = −J[σ x ⊗ σ x + h z (σ z ⊗ I + I ⊗ σ z )/2 + h x (σ x ⊗ I + I ⊗ σ x )/2] represents the local interaction between nearest neighbour spins.
Since the Hamiltonian density h keeps an overall energy scale factor J which tends to infinity in the scaling limit, we need to rescale the Trotter time step as dt = 0.05/J in order to keep small the related Trotter error. With this protocol, time in the iTEBD simulations is naturally measured in units of J. As a consequence, to reach a given iTEBD time, we need to perform a number of Trotter steps which increases linearly with J, so approaching the scaling limit more and more Trotter steps are necessary. In practice, we focused on J = 5, 10, 20; therefore, in order to reach a t max 10J, the algorithm performed ∼ 1000, ∼ 2000, and ∼ 4000 Trotter iterations, respectively.
Approaching the field theory limit, the quench in the transverse h z direction becomes smaller and smaller; however, J becomes larger and larger in such a way that the mass M = 2J|1 − h z | is kept fixed. Therefore the entanglement entropy increases with time, and the auxiliary dimension χ has to be dynamically updated in order to optimally control the truncation error. At each local step, all the Schmidt vectors corresponding to singular values larger than λ min = 10 −12 are retained. For practical reasons, to prevent the algorithm from getiting stuck in neverending computations, this condition is relaxed when χ reaches the maximal value χ max = 1024. Because of the upper bound χ max , the truncation procedure is the main source of error of the algorithm.
Let us note that for the above reasons it is extremely CPU-time consuming to reach the continuum limit of any iTEBD lattice calculation. Each simulation was executed on a GenuineIntel Core i7-4790 3.60GHz processor with 32GB of RAM. After reaching the maximum value of the auxiliary dimension χ max , the CPU-time needed for a Singular Value Decomposition (SVD) of a matrix with dimensions 2048 × 2048 is ∼ 200sec. Therefore, the total CPU-time needed to reach t max 10J is between ∼ 5 days and 18 days, depending on the value of J.