From spin chains to real-time thermal field theory using tensor networks

One of the most interesting directions in theoretical high-energy physics is understanding dynamical properties of collective states of quantum field theories. The most elementary tool in this quest are retarded equilibrium correlators governing the linear response theory. In the present letter we examine tensor networks as a way of determining them in a fully ab initio way in a class of (1+1)-dimensional quantum field theories arising as infra-red descriptions of quantum Ising chains. We show that, complemented with signal analysis using the Prony method, tensor networks calculations for intermediate times provide a powerful way to explore the structure of singularities of the correlator in the complex frequency plane.

Introduction and motivation.-Much of the progress in quantum field theory (QFT) to date has been driven by challenges posed by quantum chromodynamics (QCD)the theory of strong interactions.One such major challenge, which inspires this work, is the quest for an ab initio description of ultrarelativistic heavy-ion collisions, which probe collective states (non-zero energy density) of strong interactions in a non-equilibrium setting.The current phenomenological paradigm for these processes is based on assuming a brief far from equilibrium phase followed by a long hydrodynamic regime.Important insights about the hydrodynamic phase (such as transport coefficients) and what directly precedes it (equilibration/hydrodynamization time scales) have been gained by studying small perturbations of equilibrium in soluble models of QCD, such as the holographic (AdS/CFT) description of strongly-coupled QFTs and kinetic theory approximation to weakly-interacting QFTs.See, e.g., ref. [1] for a birds eye view on some of these questions and developments.
For an observable O, the response of an equilibrium state to a perturbation triggered by the corresponding source J is captured by the retarded thermal two-point correlator of O where ρ β is the thermal density matrix and θ(t) the Heaviside function.Such a correlator can be thought of as a change in the expectation value of δ O(t, x) due to a Dirac delta source coupled to O and localized at time 0 and position 0. For more general sources J , it is convenient to work in Fourier space (note that we distinguish functions from their Fourier transforms only through ar-guments) in which case one has The non-analytic features of G R summarize many important features of the response.For instance, for holographic QFTs described by 2-derivative gravity duals [2][3][4], the function has simple poles which give rise to exponentially decaying terms.These terms can be used to identify hydrodynamic transport coefficients as well as the characteristic time needed to transition to a hydrodynamic regime.
In addition to single poles there may be branch cuts as in free QFTs [5] and kinetic theory [6][7][8], which often give rise to power-law behaviour in time.In the presence of branch cuts there are additional subtleties that we elaborate on later.
Beyond free and strongly-coupled QFTs and kinetic theory approximations to weakly-coupled QFTs, very little is known about time dependence of collective states in QFT.In recent years, tensor networks (TN) have been revealed as a powerful tool for studying QFT on the lattice in a fully ab initio way, suitable for simulating time dependent problems or non-zero fermion density, see, e.g., refs.[9][10][11] for recent reviews of sample developments.
In this letter we explore the use of TN, in particular matrix product operator (MPO) methods, to study thermal retarded correlators in complexified frequency space in a completely ab initio way.We focus on a class of (1+1)-dimensional QFTs arising as the infra-red (IR) description of quantum Ising models (earlier studied with TN in ref. [12]).This class includes non-integrable interacting theories hosting a spectrum of non-perturbative bound states.We cross-check our MPO numerics against exact QFT predictions using the Ising model in a parameter region (vanishing longitudinal field and transverse field close to criticality) in which the IR is described by a free massive fermion QFT.We then study the quantum Ising model in the presence of a longitudinal magnetic field.We cover both the case of interacting integrable QFT arising at vanishing transverse magnetic field, and a generic non-integrable case with both the transverse and longitudinal magnetic fields.We focus on the dimension-1 fermion bilinear operator and on a homogeneous perturbation (i.e.vanishing spatial momentum p) and use the so-called Prony method, see, e.g., ref. [13], to extract the complex frequency plane structure of correlators from the numerically-obtained realtime signal.To the best of our knowledge, our results for interacting QFTs are ab initio predictions.
Ising model and IR QFTs.-The quantum Ising model is given by the Hamiltonian where J sets an overall energy scale, L denotes the total number of spins (sites), σ j x, z are Pauli matrices at position j and g / h stand for the longitudinal / transverse (magnetic) field.It is well-known that one can define complex fermionic operators b j via the Jordan-Wigner transformation [14], in terms of which When the longitudinal field vanishes (g = 0), eq. ( 3) reduces to a free fermion Hamiltonian.
We can now define two independent Majorana fermion fields as ψ(x = ja) = π/a(b † j + b j ), and ψ(x = ja) = −i π/a(b † j − b j ), where we have introduced a lattice spacing a that we take to be a = 2/J.With this prescription, in the continuum limit a → 0 the fields anticommute with each other and, on top of this, satisfy {ψ(x), ψ(y)} = { ψ(x), ψ(y)} = 2π δ(x − y).The latter conditions ensure that their two-point correlation functions in the vacuum in the CFT regime (see below) decay precisely as 1 x−y at long distances.When L → ∞, then there exists a scaling limit such that the IR (with respect to the lattice spacing a = 2 J ) description of the Hamiltonian (3) is given by the following Majorana fermions QFT: where C ≈ 0.062.The parameters M h and M g in the QFT Hamiltonian have both dimension of mass (inverse of length) and are related to the ones in the spin chain Hamiltonian (3) via M h ≡ 2J|1−h| and M g ≡ DJ |g| 8/15  with D ≈ 5.416.The scaling limit, in which the QFT description appears, then corresponds to taking M h /J → 0 while keeping the ratio M h /M g fixed, see ref. [15] for a nice discussion of the QFT limit.In addition, when temperature β −1 is included, we require β J 1 in order to only excite the IR degrees of freedom.
If M h = M g = 0, then the Hamiltonian describes a free Majorana fermion (Ising) CFT with central charge equal to 1  2 .Operators i ψψ and σ are scalar primary Hermitian operators in this CFT of dimensions ∆ = 1 and ∆ = 1  8 respectively.On the lattice, as is apparent from the relation between Hamiltonians (3) and ( 5), these operators are proportional to σ j x (with proportionality constant equal simply to −a/π) and σ j z respectively.For a class of QFTs defined as relevant deformations of the Ising CFT, i.e. by the Hamiltonian (5), there are two kinds of deformations which give rise to integrable QFT: • M h = 0, M g = 0: A massive free fermion QFT, with the fermion mass being M h .In this theory, one can find analytic expressions for the correlators.Note that in this case the same continuum theory can be represented by the spin chain in both the ferro-(h < 1) and para-magnetic (h > 1) phase.
• M h = 0, M g = 0: This is the interacting integrable E 8 QFT [16], which has a spectrum of 8 massive and stable particles -fermionic bound states.The mass of the lightest, M 1 , is given precisely by M g and the masses of the heavier particles are shown in table I.
The general form of eq. ( 5) with both i ψψ and σ turned on gives rise to an interacting non-integrable QFT which contains stable and unstable bound states [16].In this letter we study this regime at non-zero temperature numerically using MPO methods in combination with Prony analysis (see the text around eq. ( 9) for a discussion of this data analysis method).
Retarded thermal correlators in solvable cases.-For the transverse field Ising model, i.e. g = 0 in eq.(3), and in the limit of an infinite chain L → ∞ one can use the free fermion formulation to find a simple formula for the retarded thermal correlator for the transverse magnetization σ j x .At vanishing spatial momentum, which is the regime of interest in this letter, and upon normalizing σ j x to give i ψψ in the continuum limit, the correlator reads (6) where the phases φ k satisfy tan There is an intuitive understanding of the above formula in terms of two-particle (fermion) exchanges: the integral in this correlator expresses the fact that acting with σ j x excites a continuum

Branch cut
Matsubara poles Figure 1.The analytic structure of the retarded correlator of σx in the absence of longitudinal field, see eq. ( 6) derived using the free fermion mapping.The function contains branch cuts which makes the structure ambiguous.One choice of branch cut is shown in blue and another one in red.The latter representation decouples the UV and the IR which makes it more natural for field theory and in this representation single poles related to Matsubara frequencies arise.The Prony method prefers vertical branch cuts as one can see in fig. 2.
of states, which are distinguished by the relative momentum of the two fermions from which it is built, see eq. ( 4).
The analytic structure of this function in the complex frequency space is ambiguous, as noted in a similar context in ref. [7].The natural representation coming from eq. ( 6) is a branch cut stretching from the IR scale, set by twice the fermion mass M h , to the ultra-violet (UV) scale, set by the lattice spacing a.This is shown in blue in fig. 1.However, one can deform this branch cut to arrive at a representation which is more natural from a QFT perspective.In this representation, the UV and the IR are separated, and there are extra poles coming from the Fermi-Dirac distribution function.These are nothing but the fermionic Matsubara frequencies ω n = −i π β (2n + 1) for integer n and they are shown in red in fig. 1.Their contribution is that of a transient, exponentially suppressed effect.
Let us stress here that so far we did not take the QFT limit (β J → ∞ with β M h fixed) and, therefore, the presence of exponentially decaying terms in time evolution is just a feature of the spin chain in question.The transients do survive in the QFT limit and give rise to exponentially suppressed contributions known from QFT studies using CFT techniques and holography and motivating to a large degree present work, see, e.g., refs.[2][3][4].
To make this point apparent, let us note that in (1+1)dimensional CFT the retarded two-point function on a line at finite β is fixed by the conformal symmetry, see ref. [17] for a recent overview.For an operator of scaling dimension ∆ and at p = 0, the retarded correlator has  6) in the case M h = 0.The Prony analysis is done for many different time windows, characterized by the color legend.The results for each window is concatenated into a single plot.Thus, the fuzziness of a given pole gives a measure of robustness of our results.Branch cuts are represented by time-dependent poles, which results in rainbows.Prony correctly identifies both the transients as well as the UV branch cuts of fig. 1 single poles at which for ∆ = 1 coincide with single poles present in eq. ( 6) and that originate from Matsubara frequencies.Indeed, for a canonically normalized operator of dimension ∆ = 1 the retarded CFT correlator up to contact terms reads and has a sequence of poles at positions (7) with residues − 4π β .One can check that eq. ( 6) predicts the same residues for M h = 0 and in the limit β J → ∞.Therefore, for transient effects at g = 0 it is the value of residue rather than its position (which is always the same) that indicates being in the QFT regime.The residues indicate that the CFT regime is reached already for J β between 5 and 10.
Tensor networks setup and signal analysis.-Inorder to analyze the structure of singularities beyond exactly solvable cases, we evaluate numerically on the lattice the time and space dependent 2-point correlators of the form ∼ i tr (ρ β [σ x (t), σ x 0 ]) for t > 0, where ρ β is the thermal equilibrium state at inverse temperature β.This is achieved using standard tensor network techniques (see e.g.[18][19][20]).We use finite systems with open boundary conditions.We construct a matrix product operator (MPO) [21,22] approximation to the thermal state using the time evolved block decimation (TEBD) algorithm [23] to simulate imaginary time evolution.In this method, the evolution operator (in real or imaginary time) is written as a sequence of discrete time steps of width δt (Trotter step), and each of them is approximated by a sequence of two-body operations, which are sequentially applied on the MPO ansatz.The application of these gates would generally increase the bond dimension of the MPO, thus an approximation or truncation is performed to keep it bounded.The same method is later used to evolve the MPO in real time after having perturbed it with a local operator, in order to produce the thermal response functions [24][25][26].Thermal states of local Hamiltonians can be efficiently approximated by MPO [27,28], and the numerical error in this procedure is dominated by the truncation error induced by the real time evolution.Its magnitude can be estimated by comparing the numerical results obtained for different values of the maximal allowed bond dimension χ.In our calculations, χ is chosen large enough to make this error negligible compared to the remaining sources of numerical uncertainty and in captions of figures presenting our results we provide values of key parameters used in our simulations.
To reconstruct the analytic structure of the retarded correlator from real time data, we use Prony methods, see, e.g., ref. [13] for a review, to represent the function as a sum of complex exponentials where c k and ω k are complex numbers and N is either chosen by hand or automatically by the particular implementation used.Such functions satisfy a linear differential equation and Prony methods exploit that linearity in order to determine the frequencies ω k independently of the amplitudes c k .There are many such methods, going under different names such as ESPRIT, MUSIC, Padé, linear prediction (used in the context of TN in ref. [24]) etc.In this work we use the ESPRIT algorithm.One apparent drawback of the Prony method ( 9) is that by construction it only allows for poles in frequency space, not branch cuts.Nevertheless, branch cuts can still be discovered as they will be represented by a sequence of poles.We will apply the Prony method not to the full time range representing time-dependence of our correlation functions, but rather to a sequence of time windows.Complex frequencies which are stable across different time windows represent then poles and streaks represent branch cuts.To make it more apparent in our plots we use different colors to denote different time windows.The fuzziness of a given structure is a measure of its uncertainty.(Detailed descriptions for a quantitative estimation of this uncertainty and further details are provided in the supplemental material.)Results of the data analysis contain also very fuzzy features that we interpret as unphysical and arising as artifacts of various truncation errors, as well as limitations of the method itself.To test the promise of the method, we apply it to a case we understand well through other methods.
Fig. 2 shows an application of the Prony method to eq. ( 6) for M h = 0, i.e. the Ising CFT on a lattice.Comparing with exact predictions encapsulated by fig. 1 one sees that the Prony analysis identifies both branch cuts and multiple decaying poles lying further away in the complex plane.This can be contrasted with the standard Fourier transform, used in a related context in, e.g., ref. [29], which specializes in identifying poles that are close to the real axis.
As we have remarked before, the alignment of branch cuts and resulting pole structure is ambiguous, so Prony must implicitly make some choice here.In our experience, the method prefers to align branch cuts vertically, along the axis corresponding to the decay rate.Intuitively, such a choice is very efficient at late times because most contributions will be heavily suppressed.Indeed, this is very similar to the reasoning we outlined earlier when justifying efficiency of signal representation in terms of a sum of exponentially suppressed terms when all singularities are single poles.
In the next two sections we will use TN + Prony to predict the structure of singularities of correlators in interacting QFTs defined by Hamiltonian (5).
Singularities of the retarded thermal correlator: masses.-Fora free fermion CFT, the branch cut along the real axis at p = 0 arises from the exchange of pairs of fermions of zero net momentum.As we already indicated, turning on a non-zero longitudinal field leads to a confinement of fermions and non-perturbative formation of bound states (mesons).Meson masses are known exactly in the integrable E 8 QFT and otherwise were / can be determined numerically using truncated conformal space approach / DMRG [19,30,31].
The purpose of this section is to use TN+Prony to show how mesons enter the retarded correlator of i ψψ at non-zero temperature β −1 and vanishing spatial momentum.The main feature we see for the range of temperature we probe is that mesons appear in the correlator primarily as single poles corresponding to meson masses (and decay rate if unstable), see.figs.3 and 4.
Regarding the accuracy of our prediction, in table I we benchmark the results of a low temperature (basically vacuum) simulation extracted using TN+Prony with the analytic results for 7 out of 8 mesons existing in the E 8 QFT and note a rather remarkable agreement.In particular, in all the cases but the most massive meson we detect, our predictions match analytics within 1.5% and in three cases within a fraction of percent.
We then consider interacting non-integrable case starting with a low temperature (with respect to the mass of the lightest meson M 1 ).The result from TN+Prony is shown in fig. 3.As in the E 8 case, the frequencies agree with masses previously calculated in [32,33].However, some poles also develop an imaginary part, a signal that they are unstable.In particular, the ratio between the imaginary parts of the fourth and fifth mesons (see inset in fig. 3) from TN+Prony is 0.22 ± 0.04, which should be compared with the value obtained in ref. [34], 0.233.This feature comes from the absence of integrability and a presence of continuum of states starting from the threshold of twice the mass of the lightest meson M 1 , which includes states of energy in the neighbourhood of unstable mesons masses.
Apart from the mesons, there is a fuzzy structure appearing at a frequency slightly below M 1 as well as another pole with frequency slightly above than M 3 .This lower frequency arises due to open boundary condi- Ratios of masses of mesons extracted from TN+Prony compared to the analytical expectations for the integrable E8 theory, e.g.M2/M1 = 2 cos π 5 [16].Note that in the article we use Mn with n = 1, . . . to label masses of mesons for various values of Mg/M h , not only the E8 QFT case.  .Prony reconstruction of the retarded correlator of ψψ at the highest temperatures we reliably achieve in a nonintegrable ferromagnetic case.The vertical lines indicating masses is data from ref. [32].The dashed lines indicate differences in masses.As these were not seen in the vacuum result, this is a result of heating the system.In addition, there is a pole close to the origin, which has a positive imaginary part.The physical origin of this pole is unknown and we believe it is an artifact of our numerics.Simulation parameters: L = 200, χ = 170, Jδt = 0.02, Jtmax = 50, 2 nd order Trotter decomposition.
tions we use, which support excitations living close the edges of the chain.The boundary state energy is found by a DMRG calculation.The higher frequency lies at 2M 1 (dashed line in fig.3), which is the aforementioned threshold of the two particle continuum indexed by the relative momenta of the two particles.This should be associated with a branch cut, i.e. a time dependent pole, and we believe the latter is a source of the fuzziness.We should note that the results discussed so far in this section are very close to the vacuum.
Calculations are harder when the temperature is raised and we are forced to stop when we approach the regime where βM 1 ∼ 1.The clearest signature of thermal effects (understood as features of the correlator not present at very low temperatures) is the appearance of poles at locations corresponding to differences in masses, see fig. 4.Only M 12 = M 1 − M 2 is obvious, but there are signs of more.Up to our accuracy, we are not able to assess whether the mass or decay rates of mesons change with increasing temperature, but we can say that residues associated with these single poles in the correlators decrease.In practice, what we mean by this is that the corresponding coefficients c k in eq. ( 9) decay with temperature at fixed values of M h /M g and M g /J.It is to be expected that at sufficiently high temperatures the natural degrees of freedom become fermions (in rough analogy to deconfinement in QCD) and we leave this fascinating, yet numerically challenging problem for future work.
Singularities of the retarded thermal correlator: the lead-  ing transient.-Beyondmeson signatures, an intrinsically thermal feature of the linear response theory are single poles in the complex frequency plane that lead to features of the correlator decaying exponentially in time (transients).In fig. 1 one sees indications that a CFT prediction for the two least damped contributions to the retarded correlator of i ψψ is captured by the Prony method applied to the exact formula (6).Similar finding turns out to also hold when one applies TN to the same setup.
The question that we want to address now is what happens to singularities of the retarded correlator residing in free QFTs at locations specified by eq. ( 7) when we turn on interactions, i.e. for non-zero g (M g ).To shed light on this issue, we fix the transverse perturbation βM h = 0.5 and investigate the behaviour as the longitudinal perturbation is increased as βM g ≈ {0.27, 0.54, 1.08, 1.62}.The continuum limit is approached by increasing βJ = 6, 8, and 10.TN simulations of the correlator are then performed for a spin chain of size L = 100 in a time interval up to J t max = 10.
What we observe is that the single pole singularity governing the leading exponential decay not only survives in the presence of interactions, but, quite surprisingly, for the value of parameters we consider we do not observe significant systematic deviations from the free fermion QFT result 1 2π Im(β ω) = −1, see fig. 5. We refer to the supplemental material for details on the uncertainty estimation based on the Prony analyses and further remarks on MPO simulations in the (integrable) continuum limit.When increasing the interaction scale M g with respect to both β and M h , and also when increasing β J, the uncertainty in determining the value of the corresponding complex frequency grows up to 13 % in the ferromagnetic phase.
Finally, we want to address the question, whether we can identify the leading transient and meson states simultaneously through our Prony analyses.Such an example is shown in fig.6 for the largest value of the longitudinal perturbation.The decaying pole is naturally most visible at early time windows in the Prony analyses, while the first two meson states are most stable at late times.This demonstrates that the correlator signal contains proper information of the QFT regime and thermal effects, which both can be resolved by the Prony method.
Summary and outlook.-Thepresent letter should be viewed as a contribution to the study of non-equilibrium QFT using tensor networks (TN) in an innovative way.On the one hand, we perform simulations of real-time thermal correlators in spin chains in the regime of parameters where a QFT description is possible.On the other hand, although real time simulations with TN are usually reliable only for a finite time window, we show that combining them with a Prony analysis allows us to extract the analytic structure of the retarded two-point function in the complex frequency plane (as opposed to more standard Fourier transform, which probes only the vicinity of the real frequency axis).Such analytic struc-ture is known to be a quick and insightful way of summarizing the response of a system to perturbations.
The main results of our paper concern making predictions, in a fully ab initio way, about the response of thermal states to perturbatins, for a class of non-integrable interacting QFTs.In order to make our findings as robust as possible, we first developed an intuition about the QFT regime as probed by Prony in the case when our spin chain is equivalent to free fermions via Jordan-Wigner mapping, see figs. 1 and 2. We subsequently tested the accuracy of our numerics for a highly nontrivial case of the interacting integrable E 8 QFT, identifying in the response function, with very good accuracy, 7 out of 8 bound states (mesons) present in this setup, see table I.
In the most interesting case of non-integrable QFTs, we reproduced both the real and imaginary parts of several mesons states, see fig. 3, which were earlier found in refs.[16,[32][33][34].The temperature dependence of the structures in the complex frequency plane can be studied up to β −1 ∼ M 1 .We found that with increasing temperatures, poles start appearing at frequencies corresponding to differences in masses 4 and that the temperature starts affecting the residues of meson poles.Our results are stable with respect to variations in the system size, bond dimension and other numerical parameters, and they show quantitative agreement with known results for the spectrum.
Another effect triggered by the temperature are decaying (transient) contributions to the correlator , which in the free fermion case can be related to Matsubara frequencies [7] and which for holographic CFTs correspond to quasinormal modes of black holes in anti-de Sitter spacetimes [3,4].We have also extracted these numerically (see figs. 2, 5).Within our numerical precision, we find that, for the parameters we considered, these transients are not affected by interactions or integrability breaking.
Our results suggest a number of interesting problems that could be studied in future work.First of all, we have focused on the i ψψ correlator at zero momentum.But the numerical methods are equally suited to study the effect of non-vanishing momentum.It is entirely plausible that transients' momentum dependence is going to be significantly altered by interactions.Furthermore, it would be very interesting to push our calculations to larger temperatures and try to see qualitative change in the feature of the response (meson melting).In our problem, the IR limit is a massive phase, but in cases where the QFT interpolates between two CFTs in the UV and IR limits, it is to be expected that the position of transients changes, even at vanishing momentum [35], and it would be interesting to confirm it by applying our methods to a richer spin chain setup.The algorithms employed in this study do not exhaust the possibilities of TN, and it is an open question to further explore the applications of other TN methods (tensor network renormalization [36][37][38] and MERA [39,40], iMPS [41], tangent space methods [42]) in the context of time-dependence of (collective states of) (1+1)-dimensional QFTs.Finally, it would be very interesting to try to see if any of the features of our construction survive in the appropriately defined limit in which only a few (possibly coarse-grained) degrees of freedom survive -a setup amenable to quantum simulation.
We set βM h = 0.2 and approach the continuum limit by increasing βJ, using the values βJ = {2, 4, 8, 12, 16, 32}.Going to the continuum limit means that the IR length scale gets larger and we must make sure that it does not exceed the spin chain size.For a chain of L = 100 sites, the length scale of the mass compared to the chain is J/(LM h ), which may be too large for βJ > 12.This sequence corresponds to a series of increasing values of the transverse field h = {0.95,0.975, 0.9875, 0.991667, 0.99375, 0.996875}, i.e. approaching the critical point at h = 1 from the ferromagnetic phase.In Fig. 8 the result of the Prony analysis is shown for two selected values of βJ.The left column in Fig. 8 is based on the analytical result for the infinite system while the right column is based on MPO simulations for a L = 100 spin chain.From the analytical structure of the correlator, we know that there is a branch cuts stretching between the two points ±2M h /J (near the real axis) and another pair of cuts starting at 8 − 2M h stretching out to infinity.These figures demonstrate how the latter branch cut are approximated by Prony analysis through a nearly vertical line of poles in the lower imaginary plane.Furthermore, there are purely decaying thermodynamic poles on the imaginary axis, located at − 2π βJ (2n + 1) for n ∈ IN [2].In Fig. 8, these thermodynamic poles would be located at −1, −3, . . . on the rescaled imaginary axis.The first thermodynamic pole is for both values of βJ clearly visible in the analytical result as well as in the MPO simulation.For βJ = 8 also the second pole is visible in both results, while at βJ = 32 only in the analytical case.Overall, the analytic and tensor network picture agree very well.Since the values of the masses are small compared to the time windows, the branch cut between ±2M h cannot be clearly resolved.At βJ = 32, one would expect strong finite size effects for a L = 100 chain and indeed, the Prony result contains many spurious and erratic poles.As βJ increases, one is approaching the continuum limit and there are two explanations for why the results start to deviate.Either the chain is too small or the time window is too short.Since also the free fermion calculation suffers (which uses an infinite chain), it must be the time window.Figure 9 shows how the position of the thermodynamic poles converges to the expectation as we approach the continuum limit by increasing βJ.The error bars represent the lowest order with analytical result − 1 2π Re(βω) = 1 and the second pole at − 1 2π Re(βω) = 3 based on the free fermion calculation (blue, left panel) and MPO simulation (green, right panel).For this analysis, we have applied the Prony method on 75% of all discrete points in the time evolution and then successively shifted the analysis window towards later times.Based on the statistics obtained in this way, we calculated the mean value and standard deviation of all poles around the analytical positions.A selected example of such an analyses is show in Fig. 10 (left panel) for the first transient.Therein, the location of the first decaying pole on the imaginary axis is plotted for every time window.At early times, the identification is rather unstable, which is presumably related to the fact that higher order transients are not yet decayed and overlapping the signal.We thus choose a late enough time window to calculate the mean value and deviation of the pole (indicated by the red lines).The resulting uncertainty of this analyses is plotted in Fig. 9.The first thermodynamic pole can be identified in the whole temperature or mass range.For βJ ≤ 8 the pole can be resolved with an accuracy below 1%.For βJ = 32 the uncertainty increases to 5% (free fermions) or 18% (MPO).The values based on the analytical result and MPO simulation differ otherwise only marginally.
The second pole can be only identified in the temperature range βJ = 8 . . .32 (analytical case) or βJ = 8 . . .16 (MPO simulations) with large uncertainty.The reason for this behavior is the fact that at low βJ the branch cut stretching from the UV obscures the poles on the imaginary axis.On the other hand, at large βJ, the correlation length increases such that finite size and finite time effects become important in MPO simulations.Since the analytic results agrees with the MPO results, it is finite time rather than finite size that is responsible for the deviations.In addition, Fig. 11 shows the corresponding residues of the first transient in the continuum limit.The left panel is based on the Prony analyses for the critical Ising chain, while the right one is for the ferromagnetic phase with the same parameters as before in this section.The values of the residues and their uncertainties are calculated identically to the procedure for the pole location.An example of this method is shown in the right panel of Fig. 10, where the coefficients of the complex exponentials in Prony are multiplied with the inverse time dependence (because of their decay in time).The residue is then calculated as the mean value of sufficient late time windows.From Fig. 11 it is visible that the residues to approach the analytical value of the continuum limit (calculated from eq. ( 6)) with increasing values of βJ.Although the numerical resolution is not optimal, there seems to be a clear shift in the data between the critical case (βM h = 0) and the ferromagnetic phase (βM h = 0.2).This suggests that the Prony analyses is sensitive enough to capture the effect of the finite transverse mass, and importantly, shows that for large enough βJ the result is indeed a signal in the QFT regime.In more detail, the data in Fig. 11 also shows that the residue for very large values of βJ is not correctly captured in the MPO simulation, most likely because of finite size effects and short time windows.

The non-integrable case
In Fig. 12, examples of the Prony analysis of the retarded transverse correlation function are shown for two values of βM g .In the left column, the critical point is approached from the ferromagnetic phase (h < 1) and in the right column from the paramagnetic phase (h > 1).Similarly to the graph in the integrable case in Fig. 2, the UV branch cut is approximated by Prony through the nearly vertical line of poles at β/(2π) Re(ω) ≈ 13.In all examples, the first decaying thermodynamic pole on the imaginary axis is visible at 1 2π Im(βω) ≈ −1.For the largest value of the integrability breaking βM g ≈ 1.62, the corresponding uncertainty is larger in the ferromagnetic phase (lower left panel), which is visible as a a more blurred set of poles in the Prony picture.The second decaying pole is only partially visible, we therefore focus on the first pole and neglect meson frequencies for this particular study, since these would require much larger time intervals.The position of the extracted transient poles shown in Fig. 5 in the main text is generated with several parameters in the Prony method.In particular, we have chosen a cutoff value in the range 10 −6 ≤ ≤ 10 −4 (capturing how significant modes must be to be included) and the time interval, on which the Prony analysis is applied, is within the range 75 − 85% of the total time.As for the integrable case described above, the Prony analysis for the shifted time windows yields an average pole position and standard deviation.The overall uncertainty is estimated by taking the mean value of several combinations of these parameters, which then is plotted in Fig. 5.
For the same Prony parameters and with the methodology described in the free fermion case above, we calculate the corresponding residue r 1 of the first decaying pole.Fig. 13 shows the resulting uncertainties of |βr 1 | in dependence of the inverse temperature for the ferromagnetic (left panel) and paramagnetic phase (right panel).For a comparison, the analyses was also applied to the integrable case by setting βM g = 0 (black errorbars).Its corresponding analytical value in the continuum limit is shown as grey dashed line.Similarly to the pole locations, the uncertainty is increasing for larger values of βJ and βM g .For larger longitudinal perturbations, the residues seem to exhibit the tendency to decrease in the ferromagnetic and increase in the paramagnetic phase for larger values of βJ.The orange errorbars correspond to the situation when the transverse (βM h = 0.5) and longitudinal perturbation (βM g ≈ 0.54) are nearly at the same order.The corresponding data seem to be consistent with the integrable result, i.e. the longitudinal perturbation does not seem to influence the residue.

Figure 2 .
Figure 2. Prony reconstruction of the analytic structure for eq.(6) in the case M h = 0.The Prony analysis is done for many different time windows, characterized by the color legend.The results for each window is concatenated into a single plot.Thus, the fuzziness of a given pole gives a measure of robustness of our results.Branch cuts are represented by time-dependent poles, which results in rainbows.Prony correctly identifies both the transients as well as the UV branch cuts of fig.1

Figure 3 .
Figure 3. Prony reconstruction of the retarded correlator of ψψ close to the vacuum in a non-integrable ferromagnetic case.The vertical lines indicating masses are taken from ref. [32].The dashed lines indicate continuum threshold of 2M1 and a boundary state whose value is calculated with DMRG.The inset zooms in on the 4 th and the 5 th meson to show their imaginary parts.All values are consistent with QFT.Simulation parameters: L = 200, χ = 170, J δt = 0.02, J tmax = 50, 2 nd order Trotter decomposition.

Figure 4
Figure 4. Prony reconstruction of the retarded correlator of ψψ at the highest temperatures we reliably achieve in a nonintegrable ferromagnetic case.The vertical lines indicating masses is data from ref.[32].The dashed lines indicate differences in masses.As these were not seen in the vacuum result, this is a result of heating the system.In addition, there is a pole close to the origin, which has a positive imaginary part.The physical origin of this pole is unknown and we believe it is an artifact of our numerics.Simulation parameters: L = 200, χ = 170, Jδt = 0.02, Jtmax = 50, 2 nd order Trotter decomposition.

Figure 5 .
Figure 5. Extraction of the least damped transient pole for the non-integrable case using TN+Prony.The continuum limit is approached from the ferromagnetic phase for fixed βM h = 0.5 and for different values of βMg.

Figure 6 .
Figure 6.Prony reconstruction of the retarded correlator of ψψ in the non-integrable ferromagnetic case for the continuum limit.The circles represent vacuum meson masses (from ref.[32]) and the leading transient for the free fermion QFT case, respectively.The numerical results are in reasonable agreement with these data.Dashed lines indicate mass differences.Poles on the positive imaginary axis are attributed to numerical artefacts.Simulation parameters: L = 200, χ = 250, J δt = 0.005, J tmax = 50, 2 nd order Trotter decomposition.

Figure 8 .
Figure 8. Frequency analysis from Prony for the continuum limit in the integrable case.The figures show the complex values ω/J found by the Prony analyses in a particular time window as indicated by the color bar.The left column is based on the analytical result for the response function and the right column is obtained with MPO simulations.The imaginary axis is scaled by β/(2π), such that the first two thermodynamic poles are located at −1 and −3.Numerical parameters: L = 100, χ = 200, Jδt = 0.005, Jtmax = 10, 2 nd order Trotter.

Figure 9 .
Figure 9. Extraction of the purely decaying thermodynamic poles for the integrable case, based on Prony applied to analytical result (left) and MPO simulation (right).Error bars denote the calculated location including an uncertainty measure.The grey lines represent the correct result.As βJ increases, one is approaching the continuum limit and there are two explanations for why the results start to deviate.Either the chain is too small or the time window is too short.Since also the free fermion calculation suffers (which uses an infinite chain), it must be the time window.

Figure 10 .
Figure 10.Example of the Prony analysis for the MPS simulation of the retarded correlator at βJ = 8.The left plot shows the position of the first transient on the imaginary axis for every time window (with time increment step δt = 0.005).The right panel plots the corresponding coefficient, multiplied with the inverse time dependence exp[2πt/β].For time windows later than the vertical red line indicates, the average value shown by the horizontal red line was calculated.

1 • 1 •Figure 11 .
Figure 11.Extraction of the residue of the first transient in the continuum limit based on the free fermion integral (blue error bars) and MPO simulation (green error bars).The left panel is at the critical point and the right panel is calculated for finite transverse mass in the ferromagnetic phase.Grey lines represent the analytical result.

Figure 12 .
Figure 12.Frequency analysis from Prony for the continuum limit in the non-integrable case.The left column is in the ferromagnetic phase and the right column in the paramagnetic phase.The imaginary axis is scaled with βJ/(2π).The two rows correspond to two values of the integrability breaking βMg at βJ = 10.Numerical parameters: L = 100, χ = 200, Jδt = 0.005, Jtmax = 10, 2 nd order Trotter.

1 |Figure 13 .
Figure 13.Extracted residues from Pronys method for the non-integrable regime in the ferromagentic (left) and paramagnetic phase (right) in the continuum limit.The data are shown in dependence on βJ for several values of the integrability breaking βMg.Black errorbars denote the analytical value for the integrable case (βMg = 0).The curves are shown slightly displaced for graphical purposes.