Thermalization and chaos in a 1+1d QFT

We study aspects of chaos and thermodynamics at strong coupling in a scalar model using LCT numerical methods. We find that our eigenstate spectrum satisfies Wigner-Dyson statistics and that the coefficients describing eigenstates in our basis satisfy Random Matrix Theory (RMT) statistics. At weak coupling, though the bulk of states satisfy RMT statistics, we find several scar states as well. We then use these chaotic states to compute the equation of state of the model, obtaining results consistent with Conformal Field Theory (CFT) expectations at temperatures above the scale of relevant interactions. We also test the Eigenstate Thermalization Hypothesis by computing the expectation value of local operators in eigenstates, and check that their behavior is consistent with thermal CFT values at high temperatures. Finally, we compute the Spectral Form Factor (SFF), which has the expected behavior associated with the equation of state at short times and chaos at long times. We also propose a new technique for extracting the connected part of the SFF without the need of disorder averaging by using different symmetry sectors.


I. INTRODUCTION AND SYNOPSIS
The emergence of thermodynamics in closed quantum chaotic systems is a phenomenon which is still not fully understood. For example, in what sense does an initial high energy state evolving with time begin to resemble a thermal bath? The general belief is that correlation functions of local operators in such a state at sufficiently late times behave as they would at finite temperature. The Eigenstate Thermalization Hypothesis (ETH) has become the guiding framework for parameterizing such thermalization [1][2][3]. In particular, it posits that a generic high energy eigenstate of the Hamiltonian is effectively already a thermal bath as far as local observables are concerned. Thus, time dependent correlation functions in an eigenstate (or in a superposition of eigenstates in an energy band) should decay with time as expected at finite temperature, ultimately settling down to the equilibrium values predicted by the canonical ensemble. Some of the above intuition has been tested numerically in quantum many-body systems on the lattice, however there have been far fewer tests of it in the context of QFTs. The issue is finding reliable numerical methods for simulating the time evolution of excited QFT states. Hamiltonian truncation offers the possibility of such a method, as it directly computes the eigenstates and their eigenvalues. such eigenstates possesses, and how one should extract that information. Specifically, truncation (or any other calculational scheme) will always introduce errors in the states, errors which are more significant than the typical splitting between the states (which is exponentially small in the volume and energy). Thus, one should always consider observables which either involve averaging over states, or are naturally self-averaging, in order to extract physically meaningful quantities. This leads to a further question of which kind of averaging procedure is optimal for a particular observable, and which observables are naturally more self-averaging. Relatedly, on a practical level, how does one establish in the chaotic regimewhich is always non-perturbative-that the states one has computed are "healthy"?
In this paper we address these questions and study highly excited states in the context of 1+1d φ 4 theory: A schematic phase diagram expected for this theory is shown in Fig. 1. This theory is asymptotically free, so that the high temperature equation of state T m, √ λ is approximately that of a free scalar (this region is denoted 'QGP' in a generous analogy with the quark-gluon plasma). At low temperatures, the equation of state is Boltzmann suppressed as the theory is gapped, except in the Ising quantum critical fan near the critical coupling where the degrees of freedom are those of a free fermion. The theory is expected to be most strongly thermalizing and chaotic away from these free limits; we will therefore be most interested in the regime T √ λ m 1/V . The window in the figure indicates the region explored numerically in this work. Thermal equilibrium and outof-equilibrium properties of this model have been discussed in [4], we review them in appendix C.
To perform our computations, we use Lightcone Conformal Truncation (LCT), a Hamiltonian truncation method where one truncates on the scaling dimension, ∆, of primary operators in a conformal basis [5,6] (see appendix A). Our approach is distinguished from most numerical work on quantum many-body dynamics in that it is directly formulated in the continuum rather than on the lattice, and continuous translation invariance is preserved exactly in our numerics. Indeed, in LCT, the basis states describe the free massless scalar CFT. Formally, these degrees of freedom live in infinite volume. However, once we have added the relevant deformations (i.e. the mass and the quartic coupling), an effective volume emerges controlled by the truncation parameter ∆. We find that the effective volume of a state depends on its energy in a universal manner: V ∼ ∆/E, where E is the rest frame energy of the state. In fact, in practical terms, it is advantageous that the volume changes throughout the spectrum, as it effectively allows LCT to access high temperature (or energy density) states at any truncation, as we will explain below. Moreover, though there have been several Hamiltonian truncation approaches to thermalization in QFTs (see in particular Refs. [7][8][9][10][11][12] for applications of these methods to out-of-equilibrium dynamics), the regime T m, √ λ 1/V , which we explore in this work has not been studied as far as we are aware.
Having understood the above volume dependence, we can then proceed to check the "health" of our states. First, we observe that our spectrum satisfies Wigner-Dyson statistics, indicating that the eigenstates are strongly coupled and chaotic (this is true even at weak coupling λ/m 2 1). Second, using this chaotic spectrum, we compute the non-perturbative equation of state and various related observables for different choices of the coupling. These quantities exhibit Cardy behavior at high temperature, and at lower temperature agree with expectations for a CFT deformed by a relevant perturbation. Next, we use the data contained in the eigenstate wavefunctions: we verify that the coefficients describing typical eigenstates in our basis satisfy Random Matrix Theory (RMT) statistics, and use those eigenstates to calculate the expectation value of local observables-we find that those values are consistent with those of a free scalar in a thermal bath, as predicted by ETH. These tests establish the thermality of individual eigenstates. Note that this test of ETH, whereby we compare numerical expectation values to known analytic results for the free scalar thermal CFT, is possible because LCT allows us to access high temperature states.
We conclude by computing the simplest timedependent observable which captures both features of the equation of state, as well as features of chaos-namely, the Spectral Form Factor (SFF), and the related survival probability. Indeed, we find that it has a "slope-rampplateau" profile, with the slope determined by the equation of state and the ramp by RMT. We also introduce a 'symmetry-resolved' SFF, which gives access to the connected SFF without the need of disorder averaging.
Our numerical results for the model in Eq. (1) suggest that, taking symmetries properly into account, generic interacting QFTs thermalize and feature universal signatures of quantum chaos akin to quantum many-body systems on the lattice. At strong coupling, all high energy microstates are found to behave like thermal states (this can be contrasted to holographic QFTs where nonthermal or 'scar' states exist [13][14][15] despite strong thermalization and chaos). At weak coupling, we find instead a small collection of nonthermal (nonetheless nonperturbative) scar states.
Furthermore, we hope that this work is an initial step in a larger program to extend the reach of LCT methods towards the goal of understanding thermodynamics and non-equilibrium physics. In 1+1d, for any Lagrangian theory of scalars, fermions, and gauge fields, LCT methods allow one to calculate the spectrum. Previously, analysis of the spectrum was concentrated on the low lying states. However, as demonstrated in this work, one can use the high energy spectrum in order to determine the equation of state and related quantities. In certain cases, for example, when the theory is chiral, supersymmetric, or contains CP-violating phases, Hamiltonian methods can have an advantage over lattice Monte-Carlo in determining the equation of state. In addition to the spectrum, Hamiltonian methods also calculate the eigenfunctions of high energy states. Averaging over these eigenfunctions one can begin to approximate the thermal expectation values of local operators. Finally, one can attempt to use these eigenfunctions to calculate time dependent quantities, such as correlation functions in excited states. Such correlation functions should have emergent hydrodynamic behavior that can teach us about thermalization and transport. Though computationally more difficult, LCT technology is also available for calculating the above observables in 2+1d, where much less has been tested about chaos and thermalization. Hence, there are many opportunities for progress in the future.

II. THERMODYNAMICS
The basic idea of our approach is to use the nonperturbative energy eigenstates to create a thermal bath by directly performing microcanonical and macrocanonical sums. Because the truncation in LCT restricts us to a specific set of states in the theory, we want to understand when the distribution of these states produces sensible thermodynamic quantities. Ultimately, of course, we are interested in studying the time evolution of deviations around thermal ensembles, and so studying thermality is a necessary preliminary step.

Intrinsic Scan Over Volumes
One key conceptual point that we have to take into account is that the truncation effectively limits the internal size of each state-roughly speaking, our trunca-tion acts like an energy-dependent volume. This can be understood intuitively by the fact that basis states are states in the free scalar CFT which take the schematic form ∂ 1 φ∂ 2 φ · · · ∂ n φ|0 , with number of derivatives and fields bounded by the truncation ∆. One therefore expects these states (or operators) to have a spatial extent which scales with ∆ [16]. In fact, in their Fock space description, these basis states are multi-particle states with wavefunctions that are polynomials in the momenta of the individual particles, and the order of the polynomials is given by the number of derivatives of the fields in the operator description [6]. The resolution in momentum space of such a basis is inversely proportional to the order of the polynomial, again indicating that the maximum spatial extent should scale with ∆. Dimensional analysis then suggests that the volume of an eigenstate |E built out of such states scales as V ∼ ∆/E. By directly measuring eigenstates (see appendix B), we find that they indeed have a physical size V that smoothly depends on the energy of the state as where the dimensionless proportionality constant α depends on exactly how the size is defined, but has very little dependence on the state itself (for this reason, we will sometimes omit α in the following). In order to focus on the thermal properties of states, we relegate a detailed discussion of how we extract the size V of states to appendix B. An interesting feature of the volume-dependence (2) of our states is that at a fixed truncation ∆, the volume is intrinsically scanned over by the states themselves. As a consequence, one automatically sees states with large temperatures, for any value of the truncation parameter. Because the volumes are determined partially by dynamics rather than by hand, one might worry that most of the time they could end up being too small, but in fact for almost all states the volumes are larger than the inverse temperature. The reason is simply that V T ∼ √ EV is a nearly energy-independent constant, given by the total entropy S ∼ √ ∆ 1. Moreover, the fact that the entropy S or density of states e S is energy independent at fixed ∆ usefully implies that the truncation produces roughly an equal number of states per decade in energy, so that with a single diagonalization of the Hamiltonian at fixed ∆ one can study a large range in temperatures. By contrast, more traditional truncations at fixed volume produce states whose density is heavily weighted toward the highest energy scales in the truncation, since S(E) ∼ √ E. In addition to V 1/T , one must also ensure that V 1/m to study the QFT in the thermodynamic limit (the opposite regime V < 1/m, 1/ √ λ instead probes the UV CFT, with very small deformation leading to a local equilibration time larger than the volume). In terms of the energy density ε = E/V ∼ E 2 /∆, this implies that we expect to find thermal states in the window  ∆. The loose lower bound comes from the fact that Boltzmann suppression will make low temperature states thermalize inefficiently. In practice, we will find that the window of strongly thermalizing states is somewhat larger.
Knowing the size of states allows us to convert extensive quantities, such as energy and entropy in particular, into densities, by dividing by the volume V . Moreover, the energy-dependence of the volume of states must be taken into account when we extract standard thermodynamic quantities that are defined at fixed volume. In particular, T dS = dE + P dV = 1 + P ∂V ∂E dE + P ∂V ∂∆ d∆. (3)

Microcanonical Approach
One way to extract the entropy S at a given energy is to simply take the log of one over the density of states ρ(E), which is straightforward to calculate given the set of energy eigenvalues E i . More precisely, we can define Then, for each state, we obtain an entropy density s ≡ S/V and energy density ε ≡ E/V , which we plot in Fig. 2 for a specific choice of coupling λ/m 2 = 1.6. We find that, over multiple decades in ε, the entropy density exhibits Cardy behavior: We see deviations from this behavior at small and large values of ε. Deviations at small ε/m 2 are expected even in the absence of truncation, because of the physical mass gap, whereas deviations at large ε/m 2 are due to the fact that at high enough energies one simply runs out of states in the truncated spectrum. In fact, as one can see from Fig. 2, the Cardy-like growth of states extends to higher energies as the truncation parameter ∆ increases.
This microcanonical approach produces the entropy S(E, ∆) as a function of energy E and ∆, from which one can extract the temperature. By inspection of (3), using V = α∆/E. Therefore, by computing how the density of states varies with ∆ and E separately, one can extract a pressure P and temperature T for each state. However, in this paper we will not pursue this microcanonical equation of state any further. Empirically, we have found that it converges somewhat slowly with increasing ∆, and in practice a more efficient method for obtaining thermodynamic quantities is to construct canonical ensembles.

Canonical Approach
In a canonical approach, we first construct a partition function Z for our spectrum obtained at a fixed truncation ∆: where the approximation of Z as an integral holds when the spectrum is sufficiently dense. In the saddle point approximation, the integral is dominated by energies E = E * (β) wherẽ The parameterβ can be related to the usual inverse temperature β by comparing with (5): where ε = E/V . Moreover, P and E can be obtained from log Z(β, ∆) ∼ S(E, V (E, ∆)) −βE through Equations (8) and (9) are three equations that can be used to determine the three unknowns P, E and β. In particular, β =β − ∆ ∂ ∆ log Z ∂β log Z , P ε = 1 −β β , and s ε =β − log Z ∂β log Z (10) are independent of the coefficient α in the volume.
It is useful to consider the dimensionless entropy density s o (T ) ≡ s(T )/T , which was argued in [4] to be a thermodynamic C-function in 1+1d QFTs; that is, it is monotonically increasing between fixed points, and at fixed points it takes the constant value πc/3 in terms of the central charge c. In the present context, there is an O(1) unknown prefactor α in the volume of states (2), and therefore also in s o , so only the overall shape of s o (T ) is meaningful. By contrast, we can instead consider volumeinsensitive observables, for which this O(1) ambiguity factors out. Two natural such quantities are the dimensionless pressure-to-energy density ratio and the speed of sound c s . In order to extract c 2 s , we need to take derivatives with respect to ε at fixed V . Because V 2 = α∆ ε , along surfaces of fixed V we have d∆ = ∆ ε dε and (∂βε) dβ = (1 − ∆ ε ∂ ∆ ε)dε. It follows that In Fig. 3, we show s o , ε P , and c 2 s , as a function of T at strong coupling λ/m 2 = 1.6 for various truncations. The dimensionless entropy density s o still has a weak but visible dependence on ∆ at ∆ = 45, whereas the latter two observables appear to have converged fairly well even by ∆ = 30. At high temperatures T m, √ λ, these approach their CFT values ε P = c 2 s = 1, and receive corrections from the relevant deformations at lower temperatures. Notice that the speed of sound is subluminal in the temperatures of interest, consistent with causality.
Scaling collapse of (∂φ) 4 expectation value at strong and weak coupling, when plotted against E 2 /∆, as predicted by Eq. (13). Each point is obtained from a microcanonical average over 20 states.
While our numerical approach preserves exact boost invariance, it is not manifestly local and the speed of sound obtained from the spectrum is not automatically smaller than one. Observing c 2 s ≤ 1 therefore serves as a useful consistency check. At intermediate temperatures, the measured speed of sound is close to the one expected for a free massive scalar field (shown in dashed in Fig. 3). At lower temperatures T /m 2, c 2 s exceeds its value for a free massless scalar, as expected: the interaction λ = 1.6m 2 reduces the mass gap of the model, leading to less Boltzmann suppression at low temperatures. In principle, if the coupling λ is increased to its critical value where the model realizes the 1+1d Ising CFT, the speed of sound should in fact feature an upturn at low temperatures and return to the CFT value c 2 s = 1 for T → 0. While many aspects of the Ising IR CFT can be captured with LCT [6,17,18], we do not have enough resolution to observe the Cardy density of Ising states at low energies, which is responsible for this upturn.

Expectation Values
In addition to these standard thermodynamic quantities, the local operators of the theory provide a much larger set of probes of the system. In general, we expect that local operators in excited states should look approximately thermal, and in particular their expectation values should be given by a universal function of the energy density of the state. Because LCT diagonal-izes the Hamiltonian in a boosted frame, we must also account for the fact that operators with spin transform nontrivially when we boost back to the rest frame. If, for an operator O (J) with spin J, we define the expectation value where δ sets the size of a small window in energies, then taking into account the boost to the rest frame (see appendix B), locality of O (J) suggests we can write f O (E, ∆) as a function of energy density as Here,f O is defined so that we expect, up to numeric prefactors, that ε −J/2f O (ε) is a thermal expectation value depending on ∆ and E only through the combination ε. In Fig. 4, we demonstrate that this expected scaling collapse holds at strong and weak coupling for the J = −4 operator (∂ − φ) 4 . This collapse further confirms that the volume (2) has been correctly identified. Moreover, the small variance in these expectation values indicates that the ETH is satisfied for diagonal matrix elements.
One can compare these matrix elements to the theory prediction for the thermal expectation value. At high temperatures T m, √ λ these are fixed to leading order by the thermal expectation values in the free scalar CFT. (∂ − φ) 4 is a Virasoro descendant of the identity operator and therefore can acquire a thermal expectation value. Since it is dimension 4, one has (∂ − φ) 4 ∼ T 4 ∼ ε 2 . Figs. 5 and 6 show that this temperature dependence is indeed observed, in the same energy density region where where approximate Cardy growth was observed in Figs. 2 and 3. Finally, we would like to check the coefficient of the ε 2 dependence. Because the energy density is only determined up to the O(1) factor in the volume in Eq. (2), we cannot compute that numerical coefficient directly; however the same ambiguity arises in the expectation value of any other operator and therefore drops out when comparing two expectation values. Consider for example another dimension-4 operator, (∂ 2 − φ) 2 . One can show (see app. C) so that the thermal expectation value of the linear combination (∂ − φ) 4 − 5 8π (∂ 2 − φ) 2 vanishes at high temperatures (this can be understood from the fact that this linear combination is a Virasoro primary, up to a total derivative). Figs. 5 and 6 show that the expectation value of this linear combination indeed approximately vanishes at high temperature, at weak and strong coupling. A deviation at lower temperatures is expected due to the breaking of conformal invariance by the dimensionful couplings m 2 and λ. At weak coupling, we find that the deviation indeed agrees with the one predicted for a free massive scalar (dashed line in Fig. 5). Both are seen to be ∝ ε 2 and approximately equal for ε/m 2 1, as expected for the massless thermal scalar. The upper limit of the shaded region is where truncation effects become important. Bottom: a difference between the two expectation values is expected due to breaking of conformal invariance by the mass m 2 (and, to a lesser extent at weak coupling, λ). This difference agrees well with the analytic prediction for a free massive scalar (λ = 0, dashed). Despite many numerical tests of ETH in local quantum-many body systems in the literature, this is as far as we know the first quantitative test against a thermal expectation value predicted from theory.
In appendix C, we also discuss the expectation value of the operator φ 2 , which has unusual properties due to the fact that in the free massive theory it is the trace of the stress tensor. Statistics of eigenvalue spacing compared to Poisson and GOE Wigner-Dyson distributions, at strong coupling λ = 1.6m 2 after unfolding the spectrum. Only states with energy density 1 < ε/m 2 < 100 were kept. The inset shows the same data on a logarithmic scale.

III. EIGENVALUE AND EIGENVECTOR STATISTICS
An important signature of chaotic systems is that the distribution of eigenvalues and eigenstates within small energy bands should be well-described by RMT statistics. For the specific case of 1+1d φ 4 theory, whose Hamiltonian is time-reversal invariant, we specifically expect the eigenstates to be described by the Gaussian orthogonal ensemble (GOE).

Eigenvalue Spacings
Perhaps most famously, if one considers the spacing of adjacent eigenvalues in a given energy range (normalized by the average space δE in that range, i.e. "unfolded" [19]): the expectation is that if this theory is chaotic the distribution of these spacings should follow the GOE Wigner-Dyson distribution Fig. 7 shows the distribution of unfolded eigenvalue spacings for 1+1d φ 4 theory at strong coupling (λ = 1.6m 2 ) obtained with LCT at ∆ = 45. In order to focus specifically on the "healthy" thermal states satisfying Cardy growth in Fig. 3, we have restricted this distribution to eigenstates with 1 < ε/m 2 < 100. As we can see, the truncation data (blue circles) clearly follows the GOE distribution (solid gray line) rather than the Poisson distribution (dashed gray line) expected for integrable systems. One quantitative measure of how well RMT statistics describe the eigenvalue spacings is the average ratio of consecutive spacings, which was studied in [11] for the double Sine-Gordon model. One advantage of this quantity is that it does not require unfolding. The RMT prediction is r GOE ≈ 0.536 compared to the Poisson prediction r P ≈ 0.386. From our LCT data at λ = 1.6m 2 we obtain r ≈ 0.535 at ∆ = 40 and r ≈ 0.530 at ∆ = 45 (again restricting to 1 < ε/m 2 < 100). While this confirms that φ 4 theory features RMT physics at strong coupling, an obvious question remains: at which value of the coupling does the spectrum transition from Poisson to GOE (i.e. integrable to chaotic)? A convenient diagnostic for this crossover is the ratio of the distances of P (s) from the two distributions, such that η 1 corresponds to Poisson and η 1 corresponds to GOE Wigner-Dyson. Fig. 8 shows log(η) for couplings 0 ≤ λ/m 2 ≤ 0.5 (yaxis) at various levels of truncation 20 ≤ ∆ ≤ 35 (x-axis). At low ∆, it is unclear which distribution describes these weakly-coupled theories. However, as we increase the size of our basis, all of these couplings (apart from λ = 0) become better described by GOE (log η > 0). This indicates that φ 4 theory is chaotic for any nonzero coupling and deviations from GOE Wigner-Dyson at weak coupling are due to finite truncation effects.

Eigenvector Coefficients
In addition to the eigenvalues, we also expect the components of eigenvectors in a given energy band to be welldescribed by RMT. Following [20,21], we consider the overlap of energy eigenstates |E with basis states |i , In our case, the basis states are the free scalar CFT states. The GOE prediction for these coefficients is [22,23] P GOE (|c| 2 ) = 1 such that ∞ 0 dx P (x) = 1 and ∞ 0 dx xP (x) = σ. This is a particularly appealing probe of chaos as it can test for thermality of an individual eigenstate-by studying the statistics of i|E with |E fixed-in contrast to eigenvalue statistics which tests for thermality of a region of the spectrum.
In [11], it was found that the eigenstate coefficients do not match the RMT prediction (20) in multiple 1+1d QFTs, including φ 4 theory. Specifically, the tails of the distribution at small | E|i | 2 decay polynomially, rather than exponentially. While our LCT results reproduce this behavior, we find that this is a consequence of considering the overlaps of eigenstates with all basis states. If we instead restrict to basis states whose energy expectation values i|H|i are close to E, this polynomial behavior is eliminated, resulting in good agreement with RMT. Our interpretation of the need to restrict the average energy of the basis states is that the full Hamiltonian is not a truly random matrix and at large separation of energy scales the eigenstate coefficients are sensitive to theory-dependent dynamics. A more detailed discussion of this observation is presented in appendix D.
With this in mind, we can compute for each eigenstate |E the distribution of coefficients P (| E|i | 2 ) where only a fraction of basis states |i with average energy i|H|i closest to E are kept, and compare the result to the GOE distribution (20), with σ computed from P (| E|i | 2 ). We have specifically chosen to restrict to a fixed fraction (one-third) of the total number of basis states at a given ∆, so that the same number of coefficients is kept for each eigenstate; see appendix D for more details. Fig. 9 shows the resulting deviation of the eigenvector coefficients from the GOE prediction as a function of the energy density ε for a free (λ = 0), weakly-coupled (λ = 0.13m 2 ), and strongly-coupled (λ = 1.6m 2 ) theory, obtained at ∆ = 40. As expected, the free theory coefficients (green) are far from GOE, and we can clearly see multiple distinct sectors associated with states of different particle number. However, at both weak (orange) and strong (blue) coupling, we see a significant band of states which are well-described by RMT (||P −P GOE || → 0). At large values of ε, the eigenvector coefficients begin to deviate significantly from GOE due to truncation effects, as   Deviation of the coefficient statistics P (| E|i | 2 ) for individual eigenstates |E , from the GOE prediction. For each eigenstate |E , only one third of the coefficients | E|i | 2 corresponding to the basis states with average energy i|H|i closest to E were kept in the statistics. At strong coupling λ = 1.6m 2 , all eigenstates in the middle of the spectrum have close to random wavefunctions in the computational basis; at weak coupling λ = 0.13m 2 there are outliers.  their eigenvalues are near the effective UV cutoff set by ∆. At small values of ε, the coefficients also deviate from GOE, due to the fact that these are largely non-thermal states with low particle number.
One notable feature of these results is that while the majority of states at weak coupling with 1 ε/m 2 100 are well-described by GOE, we can also clearly see socalled "scar" states which are not (see [24][25][26] for reviews on many-body quantum scars). Importantly, though, these scar states have large enough energy and particle number density that they are not weakly-interacting. Fig. 10 shows that these same states violate ETH. These states are also visible in the coefficient entropy, discussed in appendix D and shown in Fig. 16. It would be interesting to study these particular states in more detail to determine whether they correspond to multi-particle states near threshold, which are expected to have a semiclassical description [27][28][29]. We note that this potential mechanism for quantum many-body scars appears to be much more robust than other known mechanisms in the literature, which typically require delicate cancellations between decay channels [30].
At strong coupling, however, we see that there are no scar states with 1 ε/m 2 300, and all coefficients in this range follow a narrow band well-described by GOE. This shows that at strong coupling all states away from the edges of the spectrum look thermal, as expected in generic interacting quantum many-body systems in the absence of disorder [31].

Spectral Form Factor
Hamiltonian truncation directly gives access to real time dynamics of quantum field theories, allowing the observation of thermalization of these systems at finite temperature. A standard probe of thermalization is the equilibrium two-point function of a local operator where we have represented the thermal expectation value microcanonically by summing over all eigenvalues E, E in a window containing N states centered at energies corresponding to the inverse temperature β. To make connection with the eigenvalue statistics in Sec. III, we consider a slightly simpler quantity that does not involve the matrix elements of an operator: the spectral form factor (SFF) [32,33]: The result, with sum computed at fixed ∆ and averaged over time, is shown in fig. 11. At early times, the SFF is dominated by its disconnected part where the absolute value is taken after performing the time averaging. This disconnected piece is well approximated by taking the density of states to be smooth. The Cardy density of states (4) together with the volume of  (22), for Z2even states (λ = 1.6m 2 ). All states with energy density 1 < ε = E 2 ∆ < 100 are kept in the microcanonical sum. The solid lines are the ramp and plateau expected from RMT for GOE (25). The initial slope is compared to the theoretical prediction ∼ 1 states (2) implies that the average density of state at fixed ∆ is constant in energy ρ(E, ∆) ≈ e √ EV ∼ e √ ∆ . One then finds leading to the observed 1/t 2 behavior of the slope of the SFF at early times in Fig. 11. At late times, the SFF is instead dominated by its connected part, and is commonly used as a test of quantum chaos: quantum chaotic systems are expected to have a late time averaged SFF that agrees with that of a random matrix. The SFF for random matrices has a distinctive ramp-plateau feature (or 'correlation hole'), which for the GOE universality class is given by [19] for t > t H .
∆ is the number of eigenstates and 2π/t H ∼ e − √ ∆ is the mean eigenvalue spacing in the microcanonical window (the time scale t H thus defined is the Heisenberg time). Fig. 11 shows excellent agreement of the SFF with this RMT prediction at late time.

Symmetry-resolved Spectral Form Factor
When global symmetries are present, RMT is only expected to describe the spectrum within superselection sectors. For example, the spectral form factor (SFF) shown in Fig. 11 agrees with the RMT prediction because only Z 2 -even states are considered, and all states have the same momentum (due to exact Lorentz invariance all momentum sectors are identical in our approach). It is interesting to further make use of the global symmetry, here Z 2 , to construct symmetry-resolved spectral form factors. Consider where Z + (t) = 1 N+ E+ e iE+t is the partition function in the Z 2 -even sector, and Z − (t) that in the Z 2 -odd sector. These are shown in Fig. 12. Both S ++ and S −− are regular SFFs, in the Z 2 -even and odd sectors respectively, and exhibit the slope-ramp-plateau profile as expected. Instead, S +− only features the slope, fixed by the equation of state: the plateau is absent because no two eigenvalues of opposite Z 2 parity are expected to match, and the ramp is absent because eigenvalues with different quantum numbers do not repel.
This observation allows for the construction of a linear combination of symmetry-resolved spectral form factors, where the slope approximately cancels: The cancelation of the slope, which is controlled by the density of states, relies on the fact that the density of states is approximately Z 2 -symmetric. Fig. 12 shows that the linear combination (27) indeed removes the slope to very high precision. This procedure therefore allows the evaluation of a form of connected SFF, without the need of disorder averaging. This is particularly appealing in theories that do not have a natural large set of couplings that one can average over; however it may also be useful for disordered theories, as disorder averages can be costly. We expect symmetry-resolved SFFs (26), (27) to be useful probes of chaos in future simulations of quantum many-body systems.
It is tempting to interpret the mild excess of the connected SFF above the GOE prediction at intermediate times 1 tm 10 ( Fig. 12 top) as possible evidence of suppression of spectral rigidity due to hydrodynamics, following [34,35] (see discussion in Sec. V). However, another possible origin is inexact cancellation between the Z 2 -even and odd density of states in (27). With the current resolution we have not been able to unambiguously establish the origin of this feature. The bottom panel in Fig. 12 shows the dependence of this feature on system size, through V = ∆/E = ∆/ε.

Survival Probability
As a final measure of real time dynamics, we will consider the 'survival probability' (28) which probes the thermalization of a basis state |i (see, e.g., [36]). While in this work we compute the survival probability by exact diagonalization of the truncated Hamiltonian, an advantage of this observable is that it can also be obtained by time evolution of the basis state through other means (e.g. [37]). The analogy with the SFF (22) is clear-the survival probability can be viewed as 'dynamical typicality' applied to the SFF [38]. Conversely, the SFF is the survival probability for the thermofield double state [39].
To reduce noise, we average (28) over several basis states with average energies i|H|i that lie in a somewhat narrow window. The result at strong coupling is shown in Fig. 13 for three different energy windows: ε/m 2 ∈ 25-35, 55-65, 115-125. The basis states in these three windows have an average particle number of 7, 9, 12, and an average scaling dimension of 35, 36, 35, respectively. The LCT results in Fig. 13 are compared to the RMT prediction (dashed lines) [19,36], which is structurally the same as for the SFF. At early times, it is dominated by the disconnected piece proportional to Eq. (24) coming from the density of states, resulting in the same 1/t 2 slope. At late times, it is dominated by the connected piece proportional to Eq. (25), leading to a similar ramp-plateau structure, though with a shallower correlation hole.

V. DISCUSSION
The results above establish that a paradigmatic nonintegrable 1+1d QFT thermalizes and displays chaos, as manifested by spectral statistics and late time dynamics following RMT. While this is the expectation for generic quantum many-body systems, most analytically tractable systems are not generic, and most numeric studies have been performed on lattices; this work suggests that continuum theories, whose microscopic structure-UV CFTs-are entirely different than a lattice, also follow this expectation. Furthermore, theoretical control over certain limits of this theory, and in particular its asymptotic freedom, allowed for quantitative tests of thermalization: the expectation value of simple operators in macroscopic high energy eigenstates agrees with thermal expectation values in the UV CFT, when the energy density of these states is large compared to the dimensionful scales of the QFT.
The Hamiltonian of a QFT of course has a rich structure far beyond that of a random matrix. In particular, locality of the underlying system should lead to the emergence of hydrodynamics at intermediate times; longwavelength modulations of hydrodynamic densities act as approximately conserved quantities that suppress spectral rigidity between eigenvalues of the Hamiltonian at energy differences larger than the Thouless energy [34]. In diffusive systems of linear size L with diffusivity D, the Thouless energy is E Th = D/L 2 -however, hydrodynamics in 1+1d QFTs is expected to be described by the KPZ universality class rather than diffusion, where instead E Th = D/L 3/2 with D entirely fixed by the equilibrium equation of state [4]. This loss of spectral rigidity can be probed through the spectral form factor [40]: at times t between the local equilibration time τ eq and the Thouless time τ Th ≡ 1/E Th , densities with wavevectors k 1/(Dt) 1/z have not fully decayed and can be used as approximate quantum numbers (z = 3/2 for KPZ and z = 2 for diffusion). In d spatial dimensions, this leads to N sectors (t) ∼ exp L d (Dt) d/z approximately decoupled sectors, and a corresponding behavior in the connected spectral form factor SFF(t) N sectors (t) t e 2S for τ eq t τ Th = L z /D [34,35]. In Fig. 11, the early time behavior of the SFF is dominated by the disconnected part, with fall-off ∼ 1/t 2 fixed by the equation of state. Removing the disconnected part is difficult without disorder averaging-the QFT under consideration (1) has a single dimensionless coupling λ/m 2 . The symmetry-resolved SFF introduced in Eq. (26) appears to efficiently capture the connected part of the SFF without averaging, and may be useful in detecting signatures of hydrodynamics, and future studies of quantum chaos in many-body systems more generally.
Emergence of hydrodynamics can also be searched for with a number of more refined observables, including equilibrium two-point functions, thermalization after a quench, entanglement growth, etc. The lower bounds on the local equilibration time τ eq of 1+1d QFTs [4] and the fact that these systems are expected to be superdiffusive further reduce the parametric window for hydrodynamic behavior; we nevertheless expect this regime to be accessible in the near future, and leave this important outstanding question for future work. Hydrodynamics in the KPZ universality class has yet to be observed numerically in non-integrable quantum many-body systems.
On the lattice, a number of methods have been developed in the past decade to study quantum dynamics without exactly diagonalizing the Hamiltonian, thereby giving access to larger systems, see Ref. [41] for a recent review. We expect that our LCT approach can be similarly improved by developing or adapting techniques specifically geared towards finite temperature dynamics. A simple yet powerful example is dynamical typicality [38], which consists in replacing the thermal state with a high (average) energy basis state. Assuming this state thermalizes, this method allows for the study of dynamics without diagonalization by solving Schrödinger equations numerically, and can be straightforwardly applied to LCT to compute e.g. the survival probability (28). The results in this work were obtained using the LCT Mathematica packages [42], which implement the calculational method described in detail in [6]. Lightcone conformal truncation is a particular Hamiltonian truncation method which uses the basis of states from a UV CFT (in this work, the free massless scalar) which can be deformed to obtain the desired QFT (1+1d φ 4 theory). The basis states are built from primary operators O in momentum space, where x ± ≡ 1 √ 2 (t ± x) and N O is an overall normalization factor. These states are formally defined in infinite volume, though, as discussed in appendix B, the resulting QFT energy eigenstates have finite volume resolution. For the free scalar, the necessary primary operators take the general form This basis is truncated to a finite-dimensional subspace by setting a maximum scaling dimension ∆ and only keeping operators with ∆ O ≤ ∆. This is equivalent to restricting to operators with at most ∆ derivatives. The QFT lightcone Hamiltonian is evaluated in this subspace, and the resulting finitedimensional matrix is numerically diagonalized to obtain the approximate QFT energy eigenstates |µ, P − . In this work, the largest basis we consider is ∆ = 45, with 44 601 Z 2 -odd states and 44 533 Z 2 -even states. The most computationally prohibitive steps in this calculation are the construction of the Hamiltonian matrix elements for the φ 4 interaction and the numerical diagonalization, which take approximately 2 hours and 6 hours on a typical laptop, respectively.

Appendix B: Size of States
In exact diagonalization studies of lattice models, the system size is a fixed external parameter given by the number of sites. In certain Hamiltonian truncation approaches to QFT, the total volume is also a fixed parameter. In LCT, the volume is instead formally infinite, but the truncation introduces an IR (as well as UV) cutoff: in particular, the highly excited states that are being used as a thermal bath in this work have a finite size, which can be measured with a simple probe operator. In this appendix, we show how this measurement produces the result quoted in Eq. (2) in the main text.
In a relativistic QFT in infinite volume, one can label states by their Lorentz-invariant momentum squared µ 2 , corresponding to the state's energy in its rest frame, and momentum P . The spectrum at any P is simply related to that at any other P by boost symmetry. To measure the volume (and shape) of a state |µ, P , one can probe it with a local operator-we will use the stress tensor, and consider its matrix elements µ, P + δ P |T νλ |µ, P with a momentum transfer δ P . The dependence of this matrix element, or form factor, on δ P indicates the shape and size of the state |µ . Since all components of the stress tensor are related by Ward identities, we will focus on a specific one below: T −− . In light-cone quantization, conserved charges are obtained by integrating over x − . One therefore has dx − T −− (x − , x + )|µ, P − = |µ, P − P − , where P − denotes the momentum of the state in the x − direction. This equation is independent of x + , by conservation. Using a Lorentz-invariant normalization of states µ , P − |µ, P − = δ µµ 2πδ(P − − P − )2P − , one finds that (B2) fixes the diagonal matrix elements of the stress tensor For finite momentum transfers q − ≡ P − − P − = 0, the matrix element is no longer fixed. In the limit where this momentum transfer is small q − P − , one can expand it as µ, P − + q − |T −− (0)|µ, P − = 2P 2 − 1 −α 2 q 2 − P 2 − + · · · (B5) A term O(q − ) is absent: it is fixed similarly to the O(1) term by using the fact that dx − T −− (x)x − is the generator of boosts. The O(q 2 − ) term is the first that is not fixed by symmetry. Its coefficient is a measure of the square  14. (a) Form factor f (q−/P−) ≡ 2 − 1 P − µ, P + q|T−−|µ, P for a state with energy density ε 10m 2 at λ = 1.6. The straight lines are fits toα 2 q 2 − /P 2 − ; this determinesα, which can be seen to increase with ∆. (b) The coefficient α ≡α/∆ is essentially independent of ∆ and the energy of the state µ, or its energy density ε ≡ µ 2 /∆. (c) The same holds in the free theory, albeit with a larger variance in the volume of states. This coefficient is measured in Fig. 14. One finds that it scales with our truncation parameter asα = ∆α, where α ∼ 1 is essentially independent of ∆, and has a very weak average dependence on the eigenstate energy µ. We therefore conclude that states with rest frame energy µ at a fixed truncation ∆ have a volume Therefore, the spectrum obtained numerically with LCT at a fixed ∆ is not at a fixed volume: instead, the volume evolves smoothly following (B7). This has implications for canonical averages over states, and must be carefully acounted for in order to compare to more conventional (fixed volume) canonical ensembles, as discussed in Sec. II. From now on, we use the notation E ≡ µ to denote the rest mass energy of a state as in the main text. Eq. (B7) motivates defining the intensive energy density of an eigenstate |E as A non-trivial check of (B7) is the scaling collapse observed in matrix elements of other local operators, when plotted against the energy density (B8). Figs. 4 and 15 show such a scaling collapse for the operators (∂φ) 4 and φ 2 , respectively. Moreover, we expect that these matrix elements should be related to thermal expectation values after compensating for the finite volume of the states, as follows. First, we boost from the frame p − = 1 in which we usually work to the rest frame p − = E √ 2 , so E|O|E rest = (E/ √ 2) −J E|O|E , where J is the spin of the operator O (J ≡ −1 for ∂ − ∼ = p − ). Then we compensate for the fact that the states are plane-wave normalized, which in the rest frame can be written as E|E ≡ (2π)2Eδ(p x − p x ). Because of the plane-wave normalization, the normalized expectation values E|O|E E|E vanish for local observables O -in physical terms, our states E have a finite size whereas the average is being taken over all of space, which is mostly vacuum. Physical operator densities can be obtained with the replacement 2πδ(p x − p x ) → V x = α∆ E : (B9) We expect this quantity to be proportional to a thermal expectation value of the density of O. Comparing with our definitions in (13), and using ε = E 2 α∆ , we see that

Appendix C: Theory Considerations
High Temperature Behavior The high temperature regime of general 1+1d QFTs was studied in [4] using conformal perturbation theory to derive strong lower bounds on the equilibration time τ eq characterizing the emergence of hydrodynamics. However, it was shown that this naive conformal perturbation theory approach breaks down for φ 4 theory, due to the lack of a thermal mass for a 1+1d massless scalar. As a result, the high temperature equation of state is non-analytic in the coupling.
For example, in free massive scalar field theory at T m the pressure is such that the dimensionless entropy density s o ≡ s T is With the addition of a φ 4 interaction, at T √ λ we instead have P = π 6 T 2 − a 0 4π 4! λT with a 0 ≈ 0.667986 [4,43]. In both cases, we see that s o is monotonically decreasing as we decrease T . At late times, the high temperature behavior of 1+1d φ 4 theory is described by the KPZ universality class, whose dissipation is entirely fixed by thermodynamics, leading to the following causality bound on the equilibration time (see [4] for more details)

Free Thermal Expectation Values
The leading high temperature expectation values for local operators are those of the free scalar CFT, which are related to vacuum correlation functions via a conformal transformation [44]. As a result, the only operators with nonzero expectation values are those built from the stress tensor (i.e. Virasoro descendants of the vacuum).
For example, the stress tensor component T −− ≡ (∂ − φ) 2 has the thermal expectation value The next operator with a nonzero thermal expectation value is the Virasoro descendant (C6) There is also a Virasoro primary V with the same scaling dimension, From the thermal expectation values we can therefore derive (note that the thermal expectation value of any overall derivative is zero) In Figs. 5 (weak coupling) and 6 (strong coupling) we confirm that these expectation values scale as ε 2 and V ≈ 0 in high energy states. At lower temperatures, the deviation from V = 0 at weak coupling agrees with the analytic expectation for the free massive scalar. We expect that the similar deviation from V = 0 at strong coupling is also due to breaking of conformality, in this case by both dimensionful couplings λ and m 2 . Comparing Figs. 5 and 6, one also finds that the expectation values O have a very weak dependence on the coupling λ at high temperatures, in qualitative agreement with expectations given that the theory is asymptotically free.
The expecation value of φ 2 is more subtle, as this is not a well-defined primary operator in the scalar CFT. However, because the Hamiltonian density H ⊃ 1 2 m 2 φ 2 , we can obtain its high temperature expectation value directly from the pressure (C10)

Appendix D: More on Eigenvector Statistics
Exact diagonalization of the Hamiltonian provides not only the spectrum of eigenstates but also their wavefunction in a computational basis {|i } -in the present approach this is the basis of primaries of the free scalar CFT. While the coefficients c iE evidently depend on the choice of this basis, one may expect eigenstates of chaotic systems to be approximately random in any reasonable computational basis. It is therefore interesting to compare statistical properties of the coefficients c iE to those of an eigenstate of a random matrix (in a fixed basis). This test of RMT is particularly appealing because it can be used to establish the thermality of an individual eigenstate. However, it is important to keep in mind that the Hamiltonian of a physical quantum many-body system has structure beyond that of a random matrix, which is reflected in its eigenvectors and eigenvalues. A simple example is the density of state, studied in Sec. II, which must be taken into account (by unfolding the spectrum) in order to observe Wigner-Dyson statistics shown in Fig. 7. More subtle deviations from RMT include the intermediate regime of hydrodynamics, discussed in Sec. V, whereby spectral rigidity is suppressed for energy windows larger than the Thouless energy. One simple probe of thermality of a given eigenstate (D1) is the coefficient entropy [20,21] S |E ≡ − i |c iE | 2 log |c iE | 2 . (D2) For chaotic physical systems, for the reasons described above, S |E is typically of the order of but smaller than its value for a GOE random matrix of the same size N , S GOE = log(0.48N ) [20,21,45]. In Fig. 16, we show the coefficient entropy for eigenstates of the 1+1d QFT in Eq. (1). At strong coupling λ = 1.6m 2 , one finds as expected that the coefficient entropy has low variance and is close but below the GOE value. At weak coupling, the variance increases and there are outliers (scar states), discussed in Sec. III. Finally, the free theory has a smaller coefficient entropy, with larger variance. We have found however that there are tests of chaos that more clearly distinguish the theory at various couplings, which we discuss next.
As a more refined test of the thermality of an eigenstate |E , one can study the entire distribution of its Distribution of eigenvector coefficients |ciE| 2 = | i|E | at strong coupling. Performing the statistics over all coefficients (orange), one finds poor agreement with the GOE distribution (gray). The apparent polynomial tail at larger values of |ciE| 2 was discussed in [11]. Restricting to eigenvectors |E with energy density ε = E/V in a small window partially suppresses this tail (yellow). If one further restricts the basis states to have average energy density εi = i|H|i /V in the same window, the tail is further suppressed and the distribution agrees well with GOE (blue). For these three sets of data one has ||P − PGOE||1 = 0.55 (orange), 0.36 (yellow), and 0.052 (blue).
coefficients P (|c iE | 2 ), and compare it to the distribution in RMT, Eq. (20). In Ref. [11], it was found in several 1+1d QFTs that the distribution of coefficients did not agree with the GOE prediction, despite strong signatures of quantum chaos in other tests (such as eigenvalue spacing). However, as discussed above, the coefficient statistics of a QFT or any local quantum many-body system are not expected to converge to those of a random matrix, due to (among other things) a nontrivial equation of state and emergence of hydrodynamics at intermediate scales. It would be interesting to derive a theoretical expectation for the distribution P (|c iE | 2 ) that takes some of these effects into account, perhaps along the lines of what has been achieved with the spectral form factor [34,35] discussed in Sec. V; or alternatively to establish a procedure to remove these effects and allow for direct comparison with RMT (similar to the unfolded eigenvalue spacing statistics (15), which is insensitive to the density of states and hydrodynamics). In the absence of such a systematic procedure, we shall proceed heuristically as follows: given that RMT is only expected to describe the system in small energy windows (smaller than the Thouless energy), one should perform statistics only over correspondingly close eigenstates and basis states. Fig. 17 shows that this intuition is largely correct. The coefficient statistics deviates significantly from the GOE expectation. However, if one restricts both the eigenstate |E and basis state |i to have (average) energy density in a fixed window, one finds remarkably good agreement with RMT. This motivates the following test of thermality of individual eigenstates: for each eigenstate |E , one studies the distribution of eigenvec- for individual eigenstates |E , from the GOE prediction. For each eigenvector |E , one performs the statistics over the coefficients | E|i | 2 , keeping only a ratio r of all basis states |i that have average energy i|H|i closest to E. Filtering basis states in this way (i.e. r < 1) yields statistics that are more consistent with GOE.
tor coefficients |c iE | 2 = | i|E | 2 keeping only the fraction 0 < r ≤ 1 of basis states |i with average energy i|H|i closest to E. Fig. 18 shows the result of this procedure, at different values of r. Especially at smaller energy densities 0.3 ε/m 2 50, one finds that eigenstates whose coefficient statistics strongly deviate from RMT in fact agree well with RMT if one restricts the basis states over which the statistics is performed. We interpret this result as indicating that a large window of energy densities (conservatively, 1 ε/m 2 100) are chaotic -a conclusion further supported by the fact that the ETH is satisfied in these states, as shown in Figs. 4 and 15. Furthermore, this test of thermality of individual eigenstates sharply distinguishes free, weakly coupled, and strongly interacting theories, as shown in Fig. 9.