## Lack of thermalization in (1+1)-d QCD at large Nc

## Axel Cort´ es Cubero ∗ † , and Neil J Robinson †

- ∗

Institute for Theoretical Physics, Center for Extreme Matter and Emergent Phenomena, Utrecht University, Princetonplein 5, 3584 CC Utrecht, The Netherlands † Institute for Theoretical Physics, University of Amsterdam, Postbus 94485, 1090 GL Amsterdam, The Netherlands

Abstract. Motivated by recent works aimed at understanding the status of equilibration and the eigenstate thermalization hypothesis in theories with confinement, we return to the 't Hooft model, the largeN c limit of (1+1)-d quantum chromodynamics. This limit has been studied extensively since its inception in the mid-1970s, with various exact results being known, such as the quark and meson propagators, the quark-antiquark interaction vertex, and the meson decay amplitude. We then argue this model is an ideal laboratory to study non-equilibrium phenomena, since it is manifestly non-integrable, yet one retains a high level of analytic control through largeN c diagrammatics. We first elucidate what are the non-equilibrium manifestations of the phenomenon of largeN c volume independence. We then find that within the confined phase, there is a class of initial states that lead to a violation of the eigenstate thermalization hypothesis, i.e. the system never thermalizes. This is due to the existence of heavy mesons with an extensive amount of energy, a phenomenon that has been numerically observed recently in the quantum Ising chain.

Keywords. Thermalization, Eigenstate thermalization hypothesis, largeN c volume independence, 't Hooft model

## 1. Introduction

Interest in non-equilibrium quantum many-body systems has undergone a resurgence in the past decade, both from the point of view of experiment [1] and theory (see, e.g., [2] and the accompanying review articles). The condensed matter physics community has been motivated by advances in the field of cold atomic gases that rendered the nonequilibrium dynamics of isolated quantum systems observable [1, 3, 4]. At the same time, the high-energy physics community was spurred on by experimental studies of the quark-gluon plasma, first at the Relativistic Heavy Ion Collider, and later at the Large Hadron Collider [5-8]. In both communities it was rather rapidly realized that a number of fundamental issues were not as well understood as one might have thought, including the questions of how, when and why isolated quantum systems thermalize. A number of seminal works addressed these questions (see [9] and references therein), elucidating the central role conservation laws play in restricting non-equilibrium dynamics. Indeed, such restrictions may be sufficient to avoid thermalization entirely [10].

Rapidly a dichotomy emerged in describing thermalization in quantum many-body systems [10,11]. On the one hand lay systems with extensively many local conservation laws, such as quantum integrable [12] and many-body localized [13-15] systems. Here the constraints imposed by the conservation laws restrict the dynamics to such an extent that thermalization is avoided [10]. Instead, after long times local expectation values in such non-equilibrium systems equilibrate and are described in terms of a generalized Gibbs ensemble [10,16,17], a generalization of the thermal ensemble to the case of many conserved quantities. Such systems retain an extensive amount of information about the initial state, thanks to the many quantities conserved under the dynamics. On the other hand are generic systems, possessing just a small (intensive) number of conserved quantities, where information about the initial state is scrambled and lost under nonequilibrium dynamics. In the long time limit, measurements in the non-equilibrium system equilibrate to those of the thermal ensemble, whose temperature is fixed by the energy of the initial state [11]. The description of non-equilibrium dynamics of generic systems is extremely challenging. Ideas that hark back to the early 1990s the eigenstate thermalization hypothesis - play a crucial role in our understanding of thermalization in isolated quantum systems [9, 18, 19]. In section 1.1, we will discuss further the eigenstate thermalization hypothesis, and the cases in which it is violated.

This integrable/generic dichotomy has been extremely successful in describing equilibration and thermalization in isolated quantum many-body systems [10,11,16,17]. However over the past couple of years it has become increasingly apparent that there are simple quantum many-body systems that fall outside this dichotomy (see, e.g., refs. [20-51]). These systems are generic, in the sense that they do not possess extensively many conserved quantities (i.e. they are non-integrable, as can be seen via level spacing statistics), but none-the-less are not compatible with the eigenstate thermalization hypothesis. As a result, these systems show anomalous non-equilibrium dynamics, such as long-lived persistent oscillations in the time-evolution of local

observables and an absence of thermalization. On the experimental side, such anomalous non-equilibrium dynamics was recently observed in a (1+1)-d lattice system of ultracold Rydberg atoms [52]. The subset of these non-integrable systems that violate the eigenstate thermalization hypothesis is discussed in detail in section 1.2.

The dynamics of thermalization, and the eigenstate thermalization hypothesis, are non-perturbative phenomena in the interaction strength. Integrable systems are known for being exactly solvable, and there are many powerful analytical tools that can be applied to the study of non-equilibrium dynamics (see, e.g., refs. [53,54]). When studying the non-equilibrium properties of generic, non-integrable systems, however, one necessarily has less analytic control and must rely on numerical methods. We are then faced with a fundamental problem: How to study strictly non-perturbative phenomena in a (necessarily) non-integrable system? Here, we will overcome this difficulty by applying tools from the largeN c expansion of low dimensional gauge theories (see ref. [55] for a review of such approaches).

## 1.1. The eigenstate thermalization hypothesis

The eigenstate thermalization hypothesis [18, 19] lies at the heart of understanding thermalization in isolated quantum many-body systems, so it is worth recapitulating some essential details. Further information, and applications of the eigenstate thermalization hypothesis, can be found in the recent review articles [56,57].

The eigenstate thermalization hypothesis aims to address the question of how and why quantum systems thermalize. Let us consider the explicit example of a pure state | Ψ 〉 time-evolving according to a local Hamiltonian H . Local observables O in the time-evolved state are given by

$$\langle O ( t ) \rangle \equiv \langle \Psi ( t ) | O | \Psi ( t ) \rangle = \sum _ { i, j } c _ { i } c _ { j } ^ { * } e ^ { i ( E _ { j } - E _ { i } ) t } \langle E _ { j } | O | E _ { i } \rangle.$$

glyph[negationslash]

Here we denote the eigenstates of the Hamiltonian with energy E i as | E i 〉 , and c i = 〈 E i | Ψ 〉 describes the projection of the state | Ψ 〉 onto the eigenstate | E i 〉 . For a generic case with no degeneracies, E i = E j ∀ i = , j glyph[negationslash]

glyph[negationslash]

$$\langle O ( t ) \rangle = \sum _ { i } | c _ { i } | ^ { 2 } \langle E _ { i } | O | E _ { i } \rangle + \sum _ { i \neq j } c _ { i } c _ { j } ^ { * } e ^ { i ( E _ { i } - E _ { j } ) t } \langle E _ { j } | O | E _ { i } \rangle.$$

The second term on the right hand side is oscillatory. Thus in the long time limit, under a stationary phase evaluation of the sum (or taking the time-average), it averages to zero and so

$$\lim _ { t \to \infty } \langle O ( t ) \rangle = \sum _ { i } | c _ { i } | ^ { 2 } \langle E _ { i } | O | E _ { i } \rangle.$$

This is known as the diagonal ensemble result for the long-time limit, as it features only diagonal matrix elements of the operator O in the eigenbasis. It would seem that eq. (3) depends sensitively on the initial state | Ψ through the superposition coefficients 〉 c i , in

the sense that small (local) changes in | Ψ 〉 can significantly change the superposition coefficients c i = 〈 E i | Ψ . Yet if the system is thermalizing, we expect eq. (2) to coincide 〉 with the microcanonical ensemble average result

$$\langle O \rangle _ { \text{MCE} } = \frac { 1 } { N _ { \Delta E } } \sum _ { \substack { i \colon \\ | E _ { i } - E _ { \Psi } | < \Delta E } } \langle E _ { i } | O | E _ { i } \rangle.$$

Here E Ψ = Ψ 〈 | H | Ψ is the energy of the state 〉 | Ψ , ∆ 〉 E is a small energy window for the microcanonical averaging, and N ∆ E is the number of states within this energy window. Why and how can this be the case?

The eigenstate thermalization hypothesis gives a set of conditions under which eq. (3) coincides with eq. (4) (when energy E Ψ is extensive). For eigenstates | E i 〉 and | E j 〉 within a given symmetry sector, matrix elements of the local observables O should satisfy (see, e.g., ref. [58])

$$\langle E _ { i } | O | E _ { j } \rangle = f _ { O } ( E ) \delta _ { i, j } + e ^ { - S ( E ) / 2 } g _ { O } ( E, \omega ) R _ { i j },$$

with E being the average energy of the eigenstates, i.e. E = ( E i + E j ) / 2, and ω being the energy difference, ω = E i -E j . Here f O ( E ) is a smooth function of the energy E for a given operator O . Off-diagonal elements are suppressed by the thermodynamic entropy S E ( ) and related to a smooth function g O ( E,ω ) of the average energy and the energy difference, as well as a random variable R ij with zero mean and unit variance.

Let us focus on the diagonal part of eq. (5), which conjectures the form of expectation values of a local operator within a given state

$$\langle E | O | E \rangle = f _ { O } ( E ).$$

That is, for a given operator O expectation values within an eigenstate are a smooth function of the energy E alone. In such a scenario, the microcanonical ensemble can be truncated to a single state (∆ E = 0 , N ∆ E = 1) in eq. (4). Such a system is then thermalizing, by construction. It is worth mentioning that expectation values in the finite volume (rather than in the thermodynamic limit) can have some finite variance from eq. (6), which should shrink to zero in the infinite volume limit. Studies that examine eq. (6) are, by necessity, usually undertaken via exact diagonalization of small systems - see the review article [57] and references therein.

1.1.1. Systems that violate the eigenstate thermalization hypothesis It so far seems to be the case that thermalizing isolated quantum many-body systems satisfy eq. (6), see [57]. There are, however, a number of examples of quantum many-body systems where the eigenstate thermalization hypothesis is violated. In the next section, section 1.2, we discuss in more detail a particular class of such systems of relevance to this work, here we will briefly mention some other examples.

The most-well studied cases of violations of the eigenstate thermalization hypothesis are those of integrable quantum systems. Following the quantum Newton's cradle

experiment [3], it was rapidly realized that such systems do not thermalize due to the presence of extensively many conservation laws [16,17,53]. In eigenstate thermalization hypothesis studies, the presence of many (hidden) conservation laws is revealed by a clear violation of eq. (6); at fixed energy E there exist eigenstates with very different expectation values of local observables, and the variance of these remains finite in the thermodynamic limit. Many-body localized systems, in which there is an emergent integrability, also violate the eigenstate thermalization hypothesis is a similar manner [13-15].

Recently, it has been realized that there exist non-integrable quantum systems that also violated the eigenstate thermalization hypothesis. Perhaps the best studied case of this so far is the PXP model, in which there are a polynomial-in-the-systemsize number of eigenstate thermalization hypothesis violating eigenstates, often called many-body quantum scars [20-26]. Other examples include certain spin models, such as generalized Affleck-Kennedy-Lieb-Tasaki chains [28, 29] and spin-1 XY models [50], as well as fracton models [30, 31] and the thin torus limit of Landau levels at certain fillings [51].

As we have already mentioned, the non-perturbative nature of these phenomena combined with the non-integrable character of these systems restricts these studies to being either purely numerical or based upon ill-controlled analytical approximations. In this work we are presented with a rare opportunity, through the use of the largeN c expansion for low dimensional gauge theories [55], to study the breaking of the eigenstate thermalization hypothesis analytically.

## 1.2. Non-equilibrium dynamics and the eigenstate thermalization hypothesis in models with confinement

There is yet another class of non-integrable models, which would generally be thought to be thermalizing [11], that display anomalous non-equilibrium behavior. Models with confinement can exhibit an absence of thermalization, unusual non-equilibrium dynamics, and a violation of the eigenstate thermalization hypothesis. Such behavior has been seen in both continuum and lattice models, with wide-ranging examples including: holographic models [32-35], certain perturbed conformal field theories [36-40], the anisotropic 2D Ising field theory [39], the Schwinger model [41-46], the extended BoseHubbard model [47], and various spin chains [36,39,48,49]. Our aim in this work will be to establish similar anomalous dynamics and an absence of thermalization in 't Hooft's model of (1+1)-d quantum chromodynamics (QCD). Before continuing to our study, we will first spend some time reminding the reader of the anomalous non-equilibrium behavior exhibited by models with confinement.

In holographic models quenches within the confined phase (i.e. those quenches where the energy injected is not sufficient to drive a confinement-deconfinement transition) have been considered in refs. [32-35]. In the holographic hard wall model, for example, on the gravity side it was found there is insufficient energy to form a

black brane (as expected to occur for thermalization) and instead long-lived oscillatory 'scattering solutions' govern the physics [32]. Similar phenomenology was observed in the AdS soliton model [33]. Later works examined a number of thermalization indicators, including those not reliant upon the formation of a black brane or black hole horizon, finding no signatures of thermalization [35]. Both local and non-local observables were observed to undergo anomalous long-lived oscillatory dynamics in these holographic confined models [35]. We do note, however, that these behaviors may be related to the strict N c = ∞ limit taken within these works, where some remnants of integrability are present. (In the N c = ∞ limit, heavy mesons have infinite lifetime in the considered holographic model, thus there are many constants of motion in this limit. Furthermore, the quench considered in Ref. [35] generates only single meson states.)

The unusual physics exhibited by holographic confined models has also been realized in a very different setting: one-dimensional quantum spin chains. The elementary excitations of such systems, spinons, can be confined by the presence of a magnetic field or long-range exchange interactions. This observation goes back to the seminal work of McCoy and Wu on the two-dimensional Ising field theory and their study of the poles of the two-particle S -matrix [59]. Two-spinon confined excitations, analogous to the mesons of QCD, have been observed in beautiful inelastic neutron scattering [60] and THz spectroscopy [61] experiments on the ferromagnet CoNb O . Identical physics 2 6 can also arise in antiferromagnets when the longitudinal magnetic field is staggered (oscillating in sign from site-to-site in the chain), as illustrated by SrCo V O 2 2 8 [62-64]. Theoretical studies of both lattice and continuum descriptions of such spin chains have observed anomalous non-equilibrium dynamics: the suppression of light-cone spreading of correlations [36, 39, 48, 49], long-lived oscillations in local observables [37-40], and an absence of thermalization [39, 40]. Importantly, these anomalies can be linked to an explicit violation of the ETH, with single meson eigenstates being distinctly nonthermal ‡ and extending far into the many-body spectrum [39,40].

## 1.3. Overview of our main results

With this discussion of the unusual non-equilibrium physics of confined models in mind, in the remainder of this work we turn our attention to (1+1)-d QCD (known as the 't Hooft model in the largeN c limit). The physical spectrum of the 't Hooft model is qualitatively similar to the confining Ising model, consisting of a tower of meson bound states (see recent eigenstate thermalization studies [39, 40] and earlier studies of the spectrum [65]). As we will show, this implies that some of the non-equilibrium aspects of the two models are very similar. Despite their similarities, the largeN c expansion gives us significantly more analytic control over the 't Hooft model (for the Ising model only numerical non-equilibrium studies have been possible [36-40]).

Taking the number of colors, N c to be large (as first introduced by 't Hooft [66]),

‡ In the sense that expectation values of local operators in single meson eigenstates do not agree with the microcanonical average (thermal) result.

a useful correspondence arises between powers of N c and number of particles in matrix elements of gauge-invariant operators. At thermal equilibrium, we use this to show that thermal corrections to expectation values of gauge-invariant operators are suppressed by higher powers of 1 /N c . This result is related to the well-known phenomenon of largeN c volume independence [67,68]. We show how this concept can be generalized to non-equilibrium physics, leading to a hierarchy of 1 /N c -suppressed terms.

We then use the 1 /N c expansion, in combination with properties of the matrix elements that appear within this expansion, to argue that quenches from finite energy density states containing few heavy, high-energy mesons do not thermalize (in the sense that the long-time limit of local observables does not approach the relevant thermal value). This statement holds even once the non-perturbative corrections that can arise in the long-time limit are taken into account. At the heart of our argument is the absence of annihilation poles in one-particle-to-one-particle matrix elements, which can be traced back to the local nature of meson excitations with respect to local observables.

We discuss our results in the context of the non-thermal states already observed in quantum magnets with confinement § . The absence of thermalization implies a violation of the eigenstate thermalization hypothesis. Our results fits into an emerging picture in which theories with confinement can exhibit unusual thermalization properties.

## 1.4. Layout

The remainder of this work proceeds as follows. In section 2.1 we remind the reader of various exact results in the 't Hooft model in the large N c limit. Section 2.2 discusses the interactions between meson excitations of the model, while in section 3 we discuss a useful particle number/powers of N c correspondence. We use this correspondence in section 4 to find a non-equilibrium generalization of the concept of largeN c volume independence, and in section 5 to discuss the non-equilibrium dynamics following a quantum quench, which leads us to describe mesons as non-thermal states in section 5.3. We conclude in section 6.

## 2. Exact results in the 't Hooft model

2.1. Planar Feynman diagrams, quark propagator and interaction vertex

We consider herein the theory of QCD in (1+1)-d described by the Lagrangian,

$$\mathcal { L } = - \frac { 1 } { 4 } G ^ { j } _ { \mu \nu \, i } G ^ { \mu \nu \, i } _ { j } + \bar { q } ^ { i } \left ( \mathrm i \gamma ^ { \mu } D _ { \mu } - m \right ) q _ { i },$$

where the indices i, j = 1 , . . . , N c run over different colors of quarks and gluons, q i is a fermionic quark field, which takes vector values in SU N ( c ) color space, and there is a SU N ( c ) matrix-valued (in the adjoint representation of SU N ( c )) gauge field A j i µ , which

§ Note that such systems can be related, via an infrared correspondence, with compactified QCD [69].

Figure 1. The Feynman rules for (7) describing (1+1)-d QCD with N c colors.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00270v2-figure-87ee1c84-c02f-4473-b837-f64d4cc4b845.jpg)

is included in the field strength and covariant derivative, respectively defined as

$$G ^ { j } _ { \mu \nu \, i } & \equiv \partial _ { \mu } A ^ { j } _ { i \, \nu } - \partial _ { \nu } A ^ { j } _ { i \, \mu } + g \, [ A _ { \mu }, A _ { \nu } ] ^ { j } _ { i } \,, \\ D _ { \mu } q _ { i } & \equiv \partial _ { \mu } q _ { i } + g A ^ { j } _ { i \, \mu } \, q _ { j }.$$

The parameters in (7) are the quark mass, m , and the strong coupling constant, g .

It is particularly convenient to study this model within the light-cone gauge, where one of the components of the gauge field is fixed to zero,

$$A _ { - } \equiv \frac { 1 } { \sqrt { 2 } } ( A _ { 0 } - A _ { 1 } ) = 0,$$

corresponding to the light-cone coordinates x ± = x ∓ = 1 √ 2 ( x 0 -x 1 ). In this gauge the Lagrangian is given by the simple form

$$\mathcal { L } = - \frac { 1 } { 2 } \text{Tr} \left ( \partial _ { - } A _ { + } \right ) ^ { 2 } - \bar { q } ( \gamma ^ { \mu } \partial _ { \mu } + m + g \gamma _ { - } A _ { + } ) q.$$

This gauge is particularly convenient because there are no Faddeev-Popov ghosts and it is easy to see that there is no vertex for gluon self interactions.

It was shown by 't Hooft that in the N c → ∞ limit, while keeping g N 2 c fixed, some aspects of this model become exactly solvable. We will refer to the N c → ∞ limit of the model defined in (7) as the 't Hooft model. Moreover, while the 't Hooft model becomes analytically tractable, the dynamics are sufficiently rich that the model is still qualitatively similar to the finite N c case. As we will see, it can be shown that the quarks in the 't Hooft model are confined, and there is an infinite tower of meson (quark-antiquark) bound states.

The largeN c limit is simplified by the fact that the leading Feynman diagrams in orders of 1 /N c have a planar structure [70]. Diagrammatically, quarks, which carry one fundamental color index, can be represented as a single line, while gluons, in the adjoint representation with two color indices, can be drawn as a double line, as pictured in figure (1). In the largeN c limit, the leading diagrams are those which can be drawn on a single plain, while the crossing of color lines is suppressed.

Figure 2. In the largeN c limit, corrections to the fully-dressed quark propagator (represented herein as a rectangular blob) in the 't Hooft model come in the form of rainbow-like diagrams, illustrated here.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00270v2-figure-2fc38afe-cc4f-4494-b966-46f89bc629b6.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00270v2-figure-7d7ed05e-d325-4541-b9d4-9a85c7d74aa3.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00270v2-figure-785a97d4-ab43-4ec7-80d9-210a5419c94c.jpg)

Figure 3. The consistency conditions for the full dressed quark propagator S p ( ), in terms of the one-particle-irreducible contributions, Σ( ), from which one obtains p eq. (12).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00270v2-figure-50f573d7-56ed-4dd5-bed4-2efa904dc11d.jpg)

In the light cone gauge, the fact that there are no gluon self-interactions [cf. (10)], combined with the simplification of considering only planar diagrams, means that the quark propagator can be computed exactly, up to 1 /N c corrections. This computation can be easily understood diagrammatically, as was shown by 't Hooft [66]. The largeN c limit implies that the only corrections to the quark propagator come from rainbow-like diagrams, involving emission and absorption of gluon lines that do not intersect each other, as shown in figure 2.

The full dressed propagator, S p ( ) can be written in terms of the bare propagator, S 0 ( p ) as

$$S ( p ) = \frac { S _ { 0 } ( p ) } { 1 + i \Sigma ( p ) S _ { 0 } ( p ) } = \frac { i p _ { - } } { 2 p _ { + } p _ { - } - m ^ { 2 } - p _ { - } \Sigma ( p ) + i \epsilon },$$

where Σ( p ) is the sum of one-particle-irreducible diagrams.

The full dressed propagator S p ( ) can be computed exactly as the solution of an integral equation that follows from the consistency condition pictured in figure 3. This

Figure 4. Bethe-Salpeter consistency condition for the quark-quark vertex (represented herein as a dark-gray circular blob). As above, the gray rectangular boxes represent full quark propagators.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00270v2-figure-fd583e9c-cb9e-4354-89e7-dd49c9932c06.jpg)

gives

$$- \mathrm i \Sigma ( p ) = \int \frac { \mathrm d k _ { - } \mathrm d k _ { + } } { ( 2 \pi ) ^ { 2 } } \frac { \mathrm i } { k ^ { 2 } } S ( p - k ) ( - 2 \mathrm i g ) ^ { 2 } N _ { c },$$

where the integral is infrared divergent, so it is necessary to introduce a low momentum cutoff λ . It can be shown [66] that the solution of (12) is

$$\Sigma ( p ) = - \frac { g ^ { 2 } N _ { c } } { \pi } \left ( \frac { \text{sgn} ( p _ { - } ) } { \lambda } - \frac { 1 } { p _ { - } } \right ),$$

and the full dressed propagator is

$$S ( p ) = \frac { - i p _ { - } } { m ^ { 2 } - g ^ { 2 } N _ { c } / \pi + 2 p _ { + } p _ { - } + g ^ { 2 } N _ { c } | p _ { - } | / \pi \lambda - i \epsilon }.$$

Due to the presence of the cutoff λ in the propagator, which has to be taken to zero, the pole in the propagator shifts to p + →∞ . This can be interpreted as the fact that the effective mass of a single quark diverges, i.e. quarks are confined and can only exist as bound states.

With knowledge of the full quark propagator at hand, it is then possible to study the full quark-quark interaction vertex, T αβ,γδ ( p, p ′ ; r ). It was shown in [71,72] that this vertex satisfies the consistency condition (pictured in figure 4),

$$T _ { \alpha \beta, \gamma \beta } ( p, p ^ { \prime }, r ) = & \frac { i g ^ { 2 } } { ( p _ { - } - p ^ { \prime } _ { - } ) ^ { 2 } } ( \gamma _ { - } ) _ { \alpha \gamma } ( \gamma _ { - } ) _ { \beta \delta } \\ & + i g ^ { 2 } N _ { c } \int \frac { d ^ { 2 } k } { ( 2 \pi ) ^ { 2 } } \frac { ( \gamma _ { - } ) _ { \alpha \sigma } ( \gamma _ { - } ) _ { \beta \lambda } } { ( k _ { - } - p _ { - } ) ^ { 2 } } S ( k ) _ { \epsilon \mu } S ( k - r ) _ { \lambda \nu } T _ { \mu \nu, \gamma \delta } ( k, p ^ { \prime } ; r ).$$

This can be simplified by introducing the function

$$\phi ( p _ { - }, p _ { - } ^ { \prime } ; r ) = \int \mathrm d p _ { + } S _ { E } ( p ) S _ { E } ( p - r ) T ( p, p ^ { \prime } ; r ),$$

where 2 γ -S E ( p ) ≡ γ -S p γ ( ) -, and T αβ γδ ; ( p, p ′ ; r ) = ( γ -) αγ ( γ -) βδ T p, p ( ′ ; r ). The relation (16) can be inverted, and the full vertex obtained in terms of φ p, p ( ′ ; r ) via

$$T ( p, p ^ { \prime } ; r ) = \frac { i g ^ { 2 } } { ( p _ { - } - p ^ { \prime } _ { - } ) ^ { 2 } } + \frac { i g ^ { 2 } N _ { c } } { \pi ^ { 2 } } \int \mathrm d k _ { - } \frac { \phi ( k _ { - }, p ^ { \prime } ; r ) } { ( k _ { - } - p _ { - } ) ^ { 2 } }.$$

The function φ p ( -, p ′ -; r ) satisfies the λ -independent Bethe-Salpeter equation (for the full derivation see [71,72]),

$$\mu ^ { 2 } \phi ( x, x ^ { \prime } ; r ) = \frac { \pi ^ { 2 } } { N _ { c } r _ { - } ( x - x ^ { \prime } ) ^ { 2 } } + \alpha \left ( \frac { 1 } { x } + \frac { 1 } { 1 - x } \right ) \phi ( x, x ^ { \prime } ; r ) + \int _ { 0 } ^ { 1 } \frac { \left [ \phi ( x, x ^ { \prime } ; r ) - \phi ( y, x ^ { \prime } ; r ) \right ] \mathrm d y } { ( x - y ) ^ { 2 } },$$

where µ 2 = r π/g N 2 2 c , α = m π/g N 2 2 c , and we switched to the more convenient variables x = p -/r -and x ′ = p ′ -/r -. Intuitively, here x can be thought of as the share of the total light cone momentum carried by the quark within the meson. The full solutions φ x, x ( ′ ; r ) can be expressed in terms of solutions of the homogeneous equation

$$\mu ^ { 2 } \phi ( x ) = \alpha \left ( \frac { 1 } { x } + \frac { 1 } { 1 - x } \right ) \phi ( x ) + \int _ { 0 } ^ { 1 } \frac { [ \phi ( x ) - \phi ( y ) ] \mathrm d y } { ( x - y ) ^ { 2 } }.$$

(Here we suppress the x ′ and r arguments.) Equation (18) has a discrete set of solutions [66], labeled by the integer k and denoted µ k , φ k ( x ). These solutions were examined numerically by 't Hooft in [66]; here we need only knowledge of some of their interesting properties. The label k corresponds to the discrete spectrum of meson particles, bound states of a quark-antiquark pair. The solutions φ k ( x ) form an orthonormal set, such that,

$$\sum _ { k } \phi _ { k } ( x ) \phi _ { k } ( x ^ { \prime } ) = \delta ( x - x ^ { \prime } ), \quad \int _ { 0 } ^ { 1 } \phi _ { n } ^ { * } ( x ) \phi _ { m } ( x ) \mathrm d x = \delta _ { n m }.$$

For large values of k , so-called heavy meson bound states, the solutions can be approximated as

$$\max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \\ \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \min \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \text{($^{2}$} )$$

The full wave function can be expressed in terms of the homogeneous solutions, φ k ( x ), as

$$\phi ( x, x ^ { \prime } ; r ) = - \sum _ { k } \frac { \pi g ^ { 2 } } { ( r ^ { 2 } - r _ { k } ^ { 2 } ) r _ { - } } \, \int _ { 0 } ^ { 1 } \text{d} y \frac { \phi _ { k } ( x ) \phi _ { k } ^ { * } ( y ) } { ( y - x ^ { \prime } ) ^ { 2 } }.$$

Combined with eq. (17), it is then possible to recover the full quark-quark scattering amplitude, T αβ,γβ ( p, p , r ′ ), see [66]. It is worth mentioning that the scattering amplitude contains 1 /λ terms, which would diverge in the λ → 0 limit. This accounts for the fact that it is impossible for a meson state to decay into free quarks, and also, as we will later see, this ensures that the scattering amplitudes between mesons are finite.

To summarize the main results quoted in this section, in the largeN c limit of QCD in (1+1)-d, one can compute exactly the quark propagator and the quark-quark interaction vertex. The quark propagator exhibits a pole at p →∞ , which is explained by the fact that the quark is confined. From the quark-quark interaction vertex one sees that only a discrete set of quark-antiquark configurations are allowed solutions, labeled by an integer k . These allowed quark-antiquark wave functions are associated with the discrete spectrum of mesonic bound states.

Figure 5. The leading contribution to the 1 → 2 meson decay process. Mesons are represented by a wiggly line, and the new effective three-meson vertex is depicted as a black square blob. The ellipsoidal dark gray blobs represent quark-quark vertices, and the gray rectangular boxes represent full quark propagators.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00270v2-figure-990523dd-c211-4362-8da1-fb91f4d0cba6.jpg)

## 2.2. Meson interactions

The spectrum of the 't Hooft model consists of an infinite tower of meson bounds states, with discrete mass values. Quarks and gluons are confined and so absent from the physical spectrum. With knowledge of the quark propagator and of quark-quark interactions, it is possible to compute the leading contributions to mesonic scattering amplitudes at large N c . The simplest interaction, first discussed in [71], is the amplitude for 1 → 2 meson decay. As shown diagrammatically in figure 5, this amplitude can be written in terms of the known quark propagator (14) and the interaction vertex (15).

The factors of λ that appear in the quark propagator, quark vertices, and arise from the integration over the loop momentum, exactly cancel to give a finite meson decay amplitude. We quote the final result for the amplitude (derived in [71]) describing the decay of a single meson (described by meson quantum number k 1 and momentum r 1 ) into two mesons (with quantum numbers k 2 and k 3 and momenta r 2 and r 3 ),

$$A _ { k _ { 1 } ; k _ { 2 }, k _ { 3 } } ( r _ { 1 } ; r _ { 2 }, r _ { 3 } ) = & \frac { 4 g ^ { 2 } \sqrt { N _ { c } } } { \sqrt { \pi } } \left [ \int _ { 0 } ^ { r _ { 2 } - } \, \mathrm d l _ { - } \, \phi _ { 1 } \left ( \frac { l _ { - } } { r _ { 1 - } } \right ) \phi _ { 2 } \left ( \frac { l _ { - } } { r _ { 2 - } } \right ) \int _ { 0 } ^ { r _ { 3 - } } \, \mathrm d p _ { - } \frac { \phi _ { 3 } \left ( \frac { p _ { - } } { r _ { 3 - } } \right ) } { [ p _ { - } \, ( l _ { - } \, - \, r _ { 2 - } ) ] ^ { 2 } } \\ & - \int _ { r _ { 2 - } } ^ { r _ { 1 - } } \, \mathrm d l _ { - } \, \phi _ { 1 } \left ( \frac { l _ { - } } { r _ { 1 - } } \right ) \phi _ { 3 } \left ( \frac { l _ { - } \, - \, r _ { 2 - } } { r _ { 3 - } } \right ) \int _ { 0 } ^ { r _ { 2 - } } \, \mathrm d p _ { - } \frac { \phi _ { 2 } \left ( \frac { p _ { - } } { r _ { 2 - } } \right ) } { ( p _ { - } \, - \, l _ { - } ) ^ { 2 } } \right ]. \ ( 2 2 )$$

Here it is important to notice that this amplitude is O (1 / √ N c ). Thus, in the strict N c →∞ limit, the mesons are completely stable and do not decay. Similarly, one can study the 2 → 2 meson scattering amplitude in terms of the exact quark propagator and interaction vertex. This is shown diagrammatically in figure 6. Its leading contribution is at order O (1 /N c ).

As we have seen, non-trivial meson scattering involving three particles happens at order 1 / √ N c , while scattering involving four particles happens at order 1 /N c . It is easy

Figure 6. A schematic depiction of the 2 → 2 meson scattering amplitude. The effective 4-meson vertex is represented by a white circular blob.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00270v2-figure-1fe13209-26a9-4258-a666-103ba3180a62.jpg)

to see a pattern: when higher numbers of particles are involved, the leading contribution to non-trivial scattering comes at higher orders of 1 /N c . A non-trivial interaction vertex involving n particles will only contribute at order 1 /N ( n -2) / 2 c . As we will explore in the next section, there is a deep correspondence between powers of 1 /N c and the number of particles involved in a process. This can also be seen at the level of matrix elements of local operators.

## 3. Particle number/power of N c correspondence

We have seen that one can diagrammatically compute the leading contributions to scattering amplitudes under a 1 /N c expansion. This approach can also be used to study matrix elements of local operators, which are necessary ingredients when studying the non-equilibrium dynamics of such systems. Physical observables have to be gauge invariant operators, for instance we can construct a meson-like operator that is given by an quark and antiquark pair connected by a path-ordered Wilson line,

$$S _ { \gamma } ( x, y ) = \bar { q } _ { j } ( x ) W _ { \gamma } ^ { j i } ( x, y ) q _ { i } ( y ),$$

where W ji γ ( x, y ) is a path-ordered exponential of the gauge field along a path, γ , which has endpoints at x and y ,

$$W _ { \gamma } ^ { j i } ( x, y ) = \mathcal { P } \exp \left [ \int _ { \gamma } \mathrm d x ^ { \mu } A _ { \mu } ^ { j i } ( x ) \right ].$$

We call operators such as S x, y ( ) 'mesonic operators', because they are gauge invariant, bilinear in the quark field q i , and create meson states when applied to the ground state.

We can also define gauge-invariant 'gluonic operators', which do not contain the quark field. These are defined in terms of the Wilson loop, given by a Wilson line along

Figure 7. Diagrammatic depiction of the one-meson matrix element.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00270v2-figure-b7b70b65-6dc5-45f3-8611-8233b3e769bf.jpg)

a closed loop C

$$W [ \mathcal { C } ] = \text{Tr} \, \mathcal { P } \exp \left [ \oint _ { \mathcal { C } } d x ^ { \mu } A _ { \mu } ^ { j i } ( x ) \right ].$$

For the sake of simplicity, in our discussion we will focus on local mesonic operators, where the quark and antiquark are at the same point,

$$S ( x ) = \bar { q } _ { i } ( x ) q ^ { i } ( x ).$$

When computing correlation functions of this operator, we will need matrix elements of the form 〈 0 | S x ( ) | n , p 1 1 ; n , p 2 2 ; . . . 〉 , which are evaluated between the ground state | 0 , 〉 and a state with a set of mesons labeled by quantum numbers n i and momenta p i .

The zero-meson matrix element (vacuum expectation value) has been computed in [73] (see also [74, 75]) and is given (in the m → 0 limit) by ‖

$$\langle 0 | S ( x ) | 0 \rangle = - \frac { N _ { c } } { \sqrt { 1 2 } } \left ( \frac { N _ { c } g ^ { 2 } } { 2 \pi } \right ) ^ { 1 / 2 } = \mathcal { O } \left ( N _ { c } \right ).$$

For finite bare quark mass 〈 0 | S x ( ) 0 | 〉 remains O ( N c ), see ref. [79]. It was shown in [71] that the matrix elements with mesons can be computed diagrammatically. The diagram corresponding to the one-meson matrix element is given in figure 7. This yields the oneparticle matrix element

$$\langle 0 | S ( x ) | n, p \rangle = \left ( \frac { N _ { c } } { \pi } \right ) ^ { 1 / 2 } \int _ { 0 } ^ { 1 } \text{d} y \, m \left ( \frac { 1 } { y } - \frac { 1 } { 1 - y } \right ) \phi _ { n } ( y ) = \mathcal { O } ( \sqrt { N _ { c } } ).$$

For the purposes of this paper, the exact details of the matrix elements with higher number of mesons are not important: we only need to point out the fact that there is a correspondence between the number of mesons in the matrix elements and the order of 1 /N c of its leading contribution. In fact, it can be shown diagrammatically that

$$\langle 0 | S ( x ) | n _ { 1 }, p _ { 1 } ; n _ { 2 }, p _ { 2 } ; \dots ; n _ { k }, p _ { k } \rangle = \mathcal { O } \left ( \frac { 1 } { N _ { c } ^ { ( k - 2 ) / 2 } } \right ).$$

‖ Note that in the light cone gauge this result is not trivial to obtain, being subtly related to boundary conditions and zero modes. Refs. [73] and [74] work in a different gauge to avoid this issue. It is also subtle in that it signals a condensate of the mesons, which breaks the (continuous) chiral symmetry. This should be forbidden by the Mermin-Wagner-Hohenberg/Coleman theorem [76-78], but the N c = ∞ limit removes the self-interactions of the Goldstone modes that would usually leads to an infrared divergence and force 〈 0 | S x ( ) 0 | 〉 = 0.

Matrix elements with mesons in both the incoming and outgoing state can be obtained from (29) by using crossing symmetry. Thus the correspondence (29) is still true, connecting the total number of particles in incoming and outgoing states to the power of 1 /N c .

One can take advantage of this correspondence, as in [71], to compute correlation functions up to some given power of 1 /N c . For example, one can compute the leading contributions to

$$M ( q ^ { 2 } ) = \int \mathrm d x \, e ^ { i q \cdot x } \langle 0 | \mathcal { T } ( S ( x ) S ( 0 ) ) | 0 \rangle$$

by inserting a complete set of meson states between the two operators (obtaining the Lehmann representation). To leading order in 1 /N c only the one-meson contribution is relevant, such that

$$M ( q ^ { 2 } ) = \sum _ { n } \frac { | \langle 0 | S ( 0 ) | n, 0 \rangle | ^ { 2 } } { q ^ { 2 } - \mu _ { n } ^ { 2 } } + \mathcal { O } \left ( \frac { 1 } { N _ { c } ^ { 2 } } \right ),$$

from which the expected behavior characteristic of asymptotic freedom can be recovered,

$$M ( q ^ { 2 } ) \sim \left ( \frac { N _ { c } } { \pi } \right ) \ln \left ( \frac { q ^ { 2 } } { \mu ^ { 2 } } \right ) \quad \text{as} \ \ q \to \infty,$$

as discussed in detail in Ref. [71].

It is important to point out that the expansion (31) is not a low-energy approximation for the correlation function, even though it only includes contributions from single mesons. As a result, one can study the ultraviolet properties of correlation functions, such as asymptotic freedom, to arbitrarily high energy scales provided the power of 1 /N c is fixed.

In the Section 5 we will see how this deep connection between the number of mesons and powers of 1 /N c has strong consequences for the thermalization properties of the 't Hooft model.

## 4. Large N c volume independence

One of the most powerful, and surprising, results in largeN c gauge theories is the phenomenon of volume independence [67,68]. Consider a gauge invariant operator O of a confining pure gauge theory (with no fermions) in d space-time dimensions, with one of the spatial dimensions compactified into a circle of circumference β . The statement of largeN c volume independence is

$$\langle O \rangle _ { \beta } = \langle O \rangle _ { \infty } \left [ 1 + \mathcal { O } \left ( \frac { 1 } { N _ { c } ^ { 2 } } \right ) \right ],$$

where 〈 . . . 〉 β represents the ground-state average of an operator in the space-time geometry R d -1 × S 1 . This statement is true as long as β is larger than any critical

value β c , where the system may enter a deconfined phase. That is, as long as the system remains in a confined phase, expectation values of gauge invariant observables remain insensitive to the compactification of one of the spatial dimensions, up to O (1 /N c ) corrections.

## 4.1. Thermal equilibrium

The existence of volume independence in a Lorentz invariant theory has immediate consequences for thermodynamics. If instead of a spatial dimension, we compactify a Euclidean time dimension to a circumference β , then that is equivalent to gaugeinvariant observables being temperature independent:

$$\langle O \rangle _ { \beta } = \frac { \text{Tr} \left [ O \, e ^ { - \beta H } \right ] } { \text{Tr} \, e ^ { - \beta H } } = \langle 0 | O | 0 \rangle \left [ 1 + \mathcal { O } \left ( \frac { 1 } { N _ { c } ^ { 2 } } \right ) \right ],$$

as long as one does not reach a deconfining critical temperature β c . We will see in the following that this statement can be modified to gauge theories including fermions, but the power of the finite-temperature corrections are modified.

The implications of thermal independence at large N c can be understood by considering the factorization properties of correlation functions of gauge-invariant operators. Let O , O 1 2 be two such operators, then their ground-state correlation functions satisfy

$$\L ^ { \bullet \, J } _ { \bullet \, J } \langle O _ { 1 } O _ { 2 } \rangle = \langle O _ { 1 } \rangle \langle O _ { 2 } \rangle \left [ 1 + \mathcal { O } \left ( \frac { 1 } { N _ { c } ^ { \alpha } } \right ) \right ],$$

where the connected piece of the correlation function is suppressed at large N c . The power, α , depends on whether the operators O 1 , O 2 are mesonic or gluonic.

We can, in fact, write down the general 1 /N c power of the leading contribution to multi-point connected correlation functions. Given a gauge theory with both quarks and gluons, we define O g i to be purely gluonic operators (such as the Wilson loop W [ C ]) and O f i to be mesonic operators (such as S x ( )). Furthermore, we impose that these gauge invariant operators are non-factorizable, i.e. they cannot be decomposed into a product of other gauge invariant operators (for instance, S 2 ( x ) is not allowed). We normalize these operators such that

$$\langle O _ { 1 } ^ { f } O _ { 2 } ^ { f } \rangle \sim N _ { c } ^ { 0 }, \quad \langle O _ { 1 } ^ { g } O _ { 2 } ^ { g } \rangle \sim \frac { 1 } { N _ { c } }.$$

It can then be shown that the leading contributions to connected correlation functions under a large N c expansion are [80],

$$\langle O _ { 1 } ^ { g } \dots O _ { n } ^ { g } O _ { 1 } ^ { f } \dots O _ { m } ^ { f } \rangle _ { \text{connected} } \sim \mathcal { O } \left ( \frac { 1 } { N _ { c } ^ { n + m / 2 - 1 } } \right ).$$

In the thermodynamic limit, the thermal expectation value (34) is equivalent to the microcanonical ensemble average, where the average is taken over states within a small

energy window ∆ E of the thermal energy E β ; recall eq. (4),

$$\langle O \rangle _ { \text{MCE} } = \frac { 1 } { N _ { \Delta E } } \sum _ { \substack { i \colon \\ | E _ { i } - E _ { \beta } | < \Delta E } } \langle E _ { i } | O | E _ { i } \rangle.$$

If β is such that the system is still in a confined phase, then all the states | E i 〉 are color-singlets. This implies that we can create all such states by acting on the ground state with a non-factorizable gauge invariant operator (or a product of them), ˆ , such E i that

$$| E _ { i } \rangle = \hat { E } _ { i } | 0 \rangle, \quad \ \langle \hat { E } _ { i } \rangle = 0.$$

We can then express the microcanonical ensemble result for the correlation function as a sum of correlation functions of the form

$$\langle \hat { E } _ { i } O \hat { E } _ { i } \rangle = \langle O \rangle \langle \hat { E } _ { i } \hat { E } _ { i } \rangle + \langle \hat { E } _ { i } O \hat { E } _ { i } \rangle _ { \text{connected} },$$

where we have used the factorization (35). The sum over index i in the first, disconnected term of eq. (40) yields the vacuum expectation value of O , which is of order N 0 c .

The connected term is suppressed as

$$\langle \hat { E } _ { i } O \hat { E } _ { i } \rangle _ { \text{connected} } \sim \frac { \langle O \rangle } { N _ { c } ^ { \alpha _ { i } } },$$

where the power, α i , depends on whether the theory contains fermionic fields or is a pure gauge theory. In a purely gluonic gauge theory both ˆ E i and O are necessarily gluonic and are, at most, non-factorizable operators such that α i = 2. This recovers the known result for largeN c volume independence [67], eq. (34). In the case where the theory contains fermions, the leading contributions will come from ˆ E i being a mesonic, non-factorizable operator and as a result α i = 1.

As long as β stays within a confining phase ¶ , all thermal corrections to expectation values of gauge-invariant operators are suppressed by powers of 1 /N c . We now want to understand if this phenomenon of volume-independence can be generalized from thermal equilibrium to non-equilibrium dynamics.

## 4.2. Generalization to non-equilibrium dynamics

A generic non-equilibrium setup involves preparing the system in an initial density matrix, ρ 0 that does not commute with the Hamiltonian of the system, H . Timeevolution is then described through the time-dependent density matrix

$$\rho ( t ) = e ^ { - i H t } \rho _ { 0 } \, e ^ { i H t },$$

¶ The separation between a deconfined and a confined phase is not straightforward in gauge theories containing massive fermions. There are different order parameters, such as the quark condensate and Polyakov loops (see for instance [81]), whose expectation values differentiates between the two phases. The critical temperature, however, is different for the each order parameter, suggesting that there is a crossover, rather than a sharp transition, between the confined and deconfined phase. For our purpose, by confined phase we mean that all the states within the energy window | E i -E β | &lt; ∆ , are color E singlets that can be created by products of non-factorizable operators.

and the evolution of physical observables is computed via

$$\langle O ( t ) \rangle = \frac { \text{Tr} \left [ O \, \rho ( t ) \right ] } { \text{Tr} \rho _ { 0 } }.$$

Here the trace can taken by summing over eigenstates of the Hamiltonian, | E i 〉 . The non-trivial time-dependence arises from the fact that the density matrix is, generally, not diagonal in this eigenstate basis. Physical observables can then be expressed as

$$\langle O ( t ) \rangle = \frac { \sum _ { i, j } \rho _ { i j } \, e ^ { - i t ( E _ { i } - E _ { j } ) } \langle E _ { j } | O | E _ { i } \rangle } { \sum _ { i } \rho _ { i i } },$$

where ρ ij = 〈 E ρ i | 0 | E j 〉 .

We now want to examine the analogue of volume independence in these nonequilibrium scenarios, where the system remains in the confined phase. This can be defined precisely by assuming that the initial conditions are such that ρ ij is non-zero only for color-singlet states, | E , i 〉 | E j 〉 . We assume the matrix ρ i,j vanishes when evaluated on colorful states above the deconfinement threshold.

Given that the sum ∑ i,j runs only over color-singlet states (which can be created by products of color-singlet, non-factorizable operators) we arrange the contributions to the time-dependent expectation value by their 1 /N c dependence glyph[negationslash]

glyph[negationslash]

$$\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \delta \end{cases}$$

The leading contribution comes from the diagonal (time-independent) term 〈 O 〉 . The rest of the terms are suppressed as

$$\langle E _ { i } | O | 0 \rangle & \sim \langle 0 | O | E _ { i } \rangle \sim \frac { \langle O \rangle } { N _ { c } ^ { \alpha _ { 1 } } }, \quad \langle E _ { j } | O | E _ { i } \rangle \sim \frac { \langle O \rangle } { N _ { c } ^ { \alpha _ { 2 } } }.$$

If we consider a purely gluonic theory, the leading contributions come from when ˆ E , E , O i ˆ j are non-factorizable gluonic operators such that α 1 = 1 and α 2 = 2. If the theory contains fermions, then the leading contributions come from mesonic, nonfactorizable operators ˆ E , E i ˆ , j such that α 1 = 1 2 and / α 2 = 1. It is then clear that non-equilibrium dynamics affects observables only at a sub-leading orders in the 1 /N c expansion, generalizing the concept of volume independence. In thermal-equilibrium, the corrections come at order 〈 O /N 〉 2 c in pure gauge theories, or at 〈 O /N 〉 c for theories with quarks. In the general non-equilibrium case, however, the highly non-diagonal matrix elements, 〈 0 | O E | i 〉 bring new lower-order corrections, at order 〈 O /N 〉 c for pure gauge theories and at 〈 O / 〉 √ N c for gauge theories with quark.

We point out that similar observations regarding the non-equilibrium manifestations of largeN c volume independence had already been made in Ref. [82], for the case of a quantum quench of the principal chiral sigma model (PCSM). The PCSM is not a

confining gauge theory, it is in fact an integrable field theory with an SU N ( c ) L × SU N ( c ) R global chiral symmetry, and with (matrix-valued) particle excitations, correspondingly labeled with left and right SU N ( c ) color indices. The arguments presented in this section are therefore not expected to apply to the PCSM. However, when one performs quantum quenches from initial states which are SU N ( c ) L × SU N ( c ) R color singlets, it was shown that contributions to the non-equilibrium dynamics are, analogously to eq. (45), suppressed by two different powers of 1 /N c . This feature, combined with the integrability of the PCSM was used in [82] to compute the full time evolution of local observables, up to the leading order in 1 /N c .

## 5. Quenches from multi-meson states

In recent works [39, 40] (see also [36-38]) it was shown that for the quantum Ising chain with both transverse and longitudinal fields, the existence of a tower of meson states, stretching far into the many-body spectrum, leads to a violation of the eigenstate thermalization hypothesis. This is due to the fact that one can create eigenstates with a finite number (e.g., one) of mesons that carry an extensive amount of energy. Then, when a physical observable is evaluated in such a states, its expectation value does not agree with the (microcanonical) thermal expectation value at the same energy density. This violation of the eigenstate thermalization hypothesis has consequences for nonequilibrium dynamics: if one performs a quantum quench in which the initial state projects strongly onto these extensive-energy few meson states, the system does not thermalize following time-evolution.

In this section we will explore a similar scenario in the 't Hooft model, by considering a quantum quench from an initial state with a finite number of mesons (i.e., zero density), but with extensive energy. The expectation values following this quench will be compared to the expected value at the corresponding finite temperature.

## 5.1. Thermal expectation value

Given the results of the previous sections, we can straightforwardly compute thermal expectation values of local operators, at leading order in 1 /N c . For the di-quark operator, (26), the thermal average is given as usual by

$$\langle S ( x ) \rangle _ { \beta } = \text{Tr} \left [ e ^ { - \beta H } S ( x ) \right ] / Z, \quad \ Z = \text{Tr} \left [ e ^ { - \beta H } \right ].$$

The trace can be computed in the basis of meson states. In this basis, the Hamiltonian can be written effectively as [70] (see also [75])

$$H = H _ { 0 } + \frac { 1 } { \sqrt { N _ { c } } } H _ { \text{int} }, \quad H _ { 0 } = \int \mathrm d p \sum _ { n } H _ { n } ^ { 0 } ( q ) \, A _ { n } ^ { \dagger } ( p ) A _ { n } ( p )$$

where A † n ( p ) is the creation operator for a meson with quantum number n , and momentum p . The meson-meson interactions come at higher orders in 1 / √ N c , as we saw in section 2.2.

To leading order in 1 /N c , the thermal expectation value is computed by considering the contributions from one-meson states only

$$\langle S ( x ) \rangle _ { \beta } = \langle 0 | S ( x ) | 0 \rangle + \sum _ { n } \int \mathrm d p W ( n, p ) \langle n, p | S ( x ) | n, p \rangle _ { \text{connected} } + \mathcal { O } \left ( \frac { 1 } { N _ { c } } \right ),$$

with the first term on the right-hand side being the vacuum expectation value, which is of order N c , as seen in eq. (27). The matrix elements involving two mesons are of order 1 (suppressed, as we expect from largeN c volume independence) and W n,p ( ) is the Bose-Einstein thermal distribution for each type of meson. Here the connected matrix element, 〈 n, p S | ( x ) | n, p 〉 connected , denotes that we have subtracted disconnected contributions, such as 〈 0 | S x ( ) 0 | 〉〈 n, p n, p | 〉 .

With an explicit expression for the thermal expectation value of S x ( ) at hand, (48), we can test if time-evolution of a given initial state with the same energy density will reproduce the thermal result, to order N 0 c , at late times.

## 5.2. Expectation value after multi-meson quench

Accordingly, we consider an initial state given by a finite set of stationary heavy mesons,

$$| \psi _ { 0 } \rangle = | n _ { 1 }, 0 ; n _ { 2 }, 0 ; \dots ; n _ { k }, 0 \rangle,$$

chosen such that the total energy E = ∑ k i =1 µ n i ∼ L is extensive and agrees with that of the thermal state for a given temperature, β . The state (49) is not an eigenstate of the full Hamiltonian, eq. (47), being instead an eigenstate of only the non-interacting (first) term. The state will therefore evolve non-trivially under | ψ t 〉 = e -i tH | ψ 0 〉 . We are interested in computing the time-dependent expectation values of local observables, such as

$$\langle S ( x ) \rangle _ { t } = \frac { \langle \psi _ { 0 } | e ^ { i t H } S ( x ) e ^ { - i t H } | \psi _ { 0 } \rangle } { \langle \psi _ { 0 } | \psi _ { 0 } \rangle }.$$

It is convenient to express this in the interaction picture,

$$\langle S ( x ) \rangle _ { t } = \frac { \langle \psi _ { 0 } | U _ { I } ^ { \dagger } ( t ) S _ { I } ( x, t ) U _ { I } ( t ) | \psi _ { 0 } \rangle } { \langle \psi _ { 0 } | \psi _ { 0 } \rangle },$$

where S I ( x, t ) = e i tH 0 S x e ( ) -i tH 0 , U I ( ) t = e i tH 0 e -i tH and H 0 is given in eq. (47). It is then possible to expand the interaction picture time-evolution operator U I ( ) t as a perturbative series in powers of 1 /N c .

The leading contribution to the perturbative expansion of (51) is obtained by taking U I ( ) t ≈ 1. This gives the time-independent expectation value

$$\langle S ( x ) \rangle _ { t } = \langle 0 | S ( x ) | 0 \rangle + \sum _ { i = 1 } ^ { k } \frac { \langle n _ { i }, 0 | S ( x ) | n _ { i }, 0 \rangle _ { c } } { \langle n _ { i }, 0 | n _ { i }, 0 \rangle } + \mathcal { O } \left ( \frac { 1 } { N _ { c } } \right ).$$

Here the first term (involving the matrix element with zero-meson states) is of order N c and the second term (with one-meson states) is order 1.

Comparing eq. (52) with eq. (48), it is immediately clear that up to order N 0 c the two do not necessarily agree. One could be tempted, therefore, to conclude that the eigenstate thermalization hypothesis is not satisfied for the initial state (49). It is, however, too early to reach this conclusion.

As we have previously discussed, the onset of thermalization is a non-perturbative phenomenon. Integrable systems reach non-thermal steady states at late times, which are at a distance O (1) from thermal values + . At first glance it would seem that the perturbative expansion of 〈 S x ( ) 〉 t , based on the expansion of the time evolution operator U I ( ), in powers of 1 t /N c , can only generate small subleading corrections to eq. (52).

We will see, however, that it is still possible (in principle) to obtain additional O (1) corrections from the perturbative expansion of U I ( ). t The only way this can happen is if some other quantity in the expansion diverges as 1 /N c is taken to 0, such that the terms in the expansion are not actually 'small', but instead O (1). We will see that this is indeed the case, as thermalization is expected to occur only at sufficiently late times, such that t α ∼ N c →∞ , for some positive power α .

5.2.1. Generating non-perturbative corrections in the long-time limit With the aim of finding O (1) corrections to eq. (52), we are interested in computing the leading contributions as t → ∞ from the 1 /N c corrections. In order for this expansion to produce O (1) terms, we would need the O (1 /N c ) correction to include a term that grows with a positive power of time, such that at late times, the expectation value would be of the form

$$\langle S ( x ) \rangle _ { t } = \langle 0 | S ( x ) | 0 \rangle + \sum _ { i = 1 } ^ { k } \frac { \langle n _ { i }, 0 | S ( x ) | n _ { i }, 0 \rangle _ { c } } { \langle n _ { i }, 0 | n _ { i }, 0 \rangle } + \frac { t ^ { \alpha _ { 1 } } } { N _ { c } } A _ { 1 } + \frac { t ^ { \alpha _ { 2 } } } { N _ { c } ^ { 2 } } A _ { 2 } + \dots, \quad ( 5 3 )$$

where α , α , . . . 1 2 are positive numbers. A simple example of how such terms that grow with time can produce O (1) contributions is, for instance, if it happens that α n = n and A n = ( -1) n Ab n n ! , in which case the expansion exponentiates at late times:

$$\langle S ( x ) \rangle _ { t } = \langle 0 | S ( x ) | 0 \rangle + \sum _ { i = 1 } ^ { k } \frac { \langle n _ { i }, 0 | S ( x ) | n _ { i }, 0 \rangle _ { c } } { \langle n _ { i }, 0 | n _ { i }, 0 \rangle } - A + A e ^ { - \frac { b } { N _ { c } } t }.$$

In this scenario of the series exponentiating at late times, the perturbative expansion generates a non-perturbative term, A , which in principle could move lim t →∞ 〈 S x ( ) 〉 t towards the thermal value 〈 S x ( ) 〉 β . We point out that a very similar scenario of exponentiation leading to non-perturbative late-time behavior has been seen in the case of a quantum quench of the Ising field theory [84] (and the corresponding lattice model [85-87]) and the sine-Gordon model [88-90].

+ See, e.g., ref. [83] for a discussion of this in the context of time-evolution in the presence of weak integrability-breaking perturbations.

5.2.2. Absence of t → ∞ non-perturbative corrections at O (1 /N c ) We will now show that it is not possible to generate the non-perturbative corrections that would lead to thermalization in the 't Hooft model. To do so, we need to understand the structure of the matrix elements involved, and how it may be possible to produce terms that grow with time.

Let us examine the time-dependence of the leading O (1 /N c ) correction to eq. (51); there are two new contributions at this order. The first comes from the O (1 /N c ) contribution to 〈 ψ 0 | S I ( x, t ) | ψ 0 〉 . This static term arises from matrix elements with states containing two particles, bringing a contribution of the form glyph[negationslash]

$$\sum _ { \substack { i, j = 1 ; \\ j \neq i } } ^ { k } \frac { \langle n _ { i }, 0 ; n _ { j }, 0 | S ( x ) | n _ { i }, 0 ; n _ { j }, 0 \rangle _ { c } } { \langle n _ { i }, 0 ; n _ { j }, 0 | n _ { i }, 0 ; n _ { j }, 0 \rangle }.$$

This is O (1 /N c ), but time-independent and thus not of the form (53) , and will not contribute at late times.

The second, more relevant, contribution arises from the O (1 /N c ) expansion of the time-evolution operator, U I ( ). t At this order in the expansion, we include the possibility for one of the mesons in the initial state to decay into a pair of mesons with opposite momenta (consistent with momentum conservation). This contribution to the expectation value is then

$$\mu _ { C } \sum _ { i, j = 1 } ^ { k } \sum _ { i ^ { \prime }, i ^ { \prime \prime } = 1 } ^ { \infty } \sum _ { j ^ { \prime }, j ^ { \prime \prime } = 1 } ^ { \infty } \int d p _ { i } d p _ { j } & \mathcal { A } _ { j ; j ^ { \prime }, j ^ { \prime \prime } } ^ { * } ( p _ { j } ) \mathcal { A } _ { i ; i ^ { \prime }, i ^ { \prime \prime } } ( p _ { i } ) \, e ^ { i t ( E _ { i ^ { \prime }, p _ { i } } + E _ { i ^ { \prime \prime }, - p _ { i } } - E _ { j ^ { \prime }, p _ { j } } - E _ { j ^ { \prime \prime }, - p _ { j } } ) } \\ & \times \frac { \langle n _ { j ^ { \prime } }, p _ { j } | S ( x ) | n _ { i ^ { \prime } }, p _ { i } \rangle _ { c } } { \langle n _ { j }, 0 | n _ { i }, 0 \rangle } \langle n _ { j ^ { \prime \prime } }, - p _ { j } | n _ { i ^ { \prime \prime } }, - p _ { i } \rangle. \\ \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \end{array}$$

Here 1 √ N c A i i ,i ; ′ ′′ ( p i ) is the amplitude for a meson with zero momentum and quantum number n i to decay into a pair of mesons with quantum numbers n i ′ and n i ′′ and momenta p i and -p i respectively. This decay process, as discussed above, is suppressed by a factor of 1 / √ N c , which we have explicitly factored out for clarity.

The expression (56) is simplified by using the factor 〈 n j ′′ , -p j | n i ′′ , - 〉 p i to eliminate one of the momentum integrals and a summation over quantum numbers, yielding

$$\frac { 1 } { N _ { c } } \sum _ { i, j = 1 } ^ { k } \sum _ { i ^ { \prime }, i ^ { \prime \prime }, j ^ { \prime } = 1 } ^ { \infty } \int \mathrm d p _ { i } & \mathcal { A } _ { j ; j ^ { \prime }, i ^ { \prime \prime } } ^ { * } ( p _ { i } ) \, \mathcal { A } _ { i ; i ^ { \prime }, i ^ { \prime \prime } } ( p _ { i } ) \, e ^ { i t ( E _ { i ^ { \prime }, p _ { i } } + E _ { i ^ { \prime \prime }, - p _ { i } } - E _ { j ^ { \prime }, p _ { i } } - E _ { i ^ { \prime \prime }, - p _ { i } } ) } \\ & \times \frac { \langle n _ { j ^ { \prime } }, p _ { i } | S ( x ) | n _ { i ^ { \prime } }, p _ { i } \rangle _ { c } } { \langle n _ { j }, 0 | n _ { i }, 0 \rangle }. \\ \dots \quad. & \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot$$

The only time dependence present in (57) is in the oscillatory phase,

$$\exp \left [ \mathrm i t \left ( E _ { i ^ { \prime }, p _ { i } } + E _ { i ^ { \prime \prime }, - p _ { i } } - E _ { j ^ { \prime }, p _ { i } } - E _ { i ^ { \prime \prime }, - p _ { i } } \right ) \right ].$$

This appears within an integral over momentum p i and summations over quantum numbers. At late times, this generically becomes a rapidly oscillating phase and, since

we are integrating over p i and summing over the quantum numbers i , i ′ ′′ , j ′ , by stationary phase arguments we find this contribution generally decays at late times as some power law in t .

The only way to obtain a behavior other than power-law decay is if there is some singularity in the integrand, specifically in the matrix elements 〈 n j 1 , p i | S x ( ) | n i 1 , p i 〉 , which could give additional non-decaying, non-trivial contributions besides the stationary phase expectation. This has indeed been observed in similar calculations of the late-time dynamics in the Ising field theory and the sine-Gordon model, see [84-90]. There singularities in the matrix elements lead to non-perturbative late-time dynamics.

At this point we note that there are significant differences between matrix elements of the Ising field theory/the sine-Gordon model (which contain singularities that can produce growing terms at large times) and those that we consider here. In particular the matrix element 〈 n j 1 , p i | S x ( ) | n i 1 , p i 〉 c does not contain any singularities that would produce any non-trivial terms. Indeed, it is well understood in the literature of integrable models in (1+1)-d where these singularities in the matrix elements come from. These are known as 'annihilation poles' which occur when a particle in the outgoing state can annihilate with an identical particle in the incoming state, which leads generally to there being a simple pole when the momentum of the particles approach each other [91,92].

It is well known that annihilation poles are always present in matrix elements which contain two-or-more particles in both the incoming and outgoing states. On the other hand, matrix elements with only a single incoming and outgoing particle do not generally have such annihilation poles. Indeed, such matrix elements only contain annihilation poles in the case where the particles are non-local with respect to the operator in question [93, 94], S x ( ) in the case at hand. That is, when particles are topological kinks or soliton-like. This is indeed the case for the particles in the Ising field theory and the sine-Gordon model, which are solitons that interpolate between different vacua and are thus non-local relative to the field operators. In the case at hand, mesonic particles are not soliton-like particles and are, in this respect, local relative to the operator S x ( ) (and other physical local operators that may be considered). Thus the matrix element 〈 n, p S | ( x ) | n, p 〉 c should not have any singularities. This is also consistent with the absence of annihilation poles in the two meson electromagnetic form factor, see [72,95] ∗ .

We can now rule out O (1 /N c ) corrections that grow at late times in the 't Hooft model. This rules out the scenario for thermalization proposed in eq. (54).

∗ Much like the 't Hooft model, the attractive regime of the sine-Gordon model contains both nonlocal particles (solitons) and local ones (breathers, neutral soliton-antisoliton pairs), relative to the fundamental bosonic field in the Lagrangian, or local functions thereof, such as exponential vertex operators. When considering a local operator, such as the vertex operator, two particle form factors display annihilation poles when both particles are solitons, but not when they are breathers. This is a direct manifestation of the effect of locality, as discussed in the text: the solitons are non-local kinks with respect to the vertex operator, while the breathers are local. It is also worth noting that the vertex operator has a non-zero vacuum expectation value, much like our local observable (27). See refs. [89, 96] for further discussions.

5.2.3. Higher order corrections Of course, one can still argue that at some higher order in 1 /N c that a 'non-perturbative' correction of the form (54) will arise. With an understanding of the mechanism under which such terms grow at late times, we can schematically understand these higher order corrections. We need only to count the number of integrals over momenta and the annihilation poles of the matrix element. The leading contributions at late times will be those arising from matrix elements with the highest number of annihilation poles.

Let us consider some examples. Working at order 1 /N 2 c , the late time leading correction would contain connected matrix elements with two particles in both the incoming and outgoing states. These matrix element contain a single annihilation pole that can contribute to growth at late times. At order 1 /N 3 c the leading term contains three particles in both incoming and outgoing states, and two annihilation poles arise. And so forth. As has already been argued in the Ising field theory and the sine-Gordon model [84-90], each new annihilation poles contributes a term that grows with higher powers at late times. It is then natural to expect contributions to (52) at late times are of the form

$$\langle S ( x ) \rangle _ { t } = \langle 0 | S ( x ) | 0 \rangle + \sum _ { i = 1 } ^ { k } \frac { \langle n _ { i }, 0 | S ( x ) | n _ { i }, 0 \rangle _ { c } } { \langle n _ { i }, 0 | n _ { i }, 0 \rangle } + \frac { 1 } { N _ { c } } A _ { 1 } + \frac { t } { N _ { c } ^ { 2 } } A _ { 2 } + \frac { t ^ { 2 } } { N _ { c } ^ { 3 } } A _ { 3 } + \dots \quad ( 5 8 )$$

Thus it is conceivable that the growing terms re-exponentiate at late times, as discussed in (54). However, as the O (1 /N c ) term does not grow in time, the growth in time is shifted by one power of 1 /N c in the 't Hooft model (as compared to Ising and sineGordon) and the non-perturbative correction arising from the exponentiation will affect only the O (1 /N c ) term.

We see, therefore, that the O (1) difference between the long-time expectation value and the thermal result cannot be overcome by the non-perturbative corrections, which modify only the 1 /N c correction, even at late times. The expression for the expectation value (52) therefore holds true, even in the t → ∞ limit. We can thus say that the thermal expectation value, (48), is generally not recovered at late times after the quantum quench. The 't Hooft model in the largeN c expansion therefore does not thermalize.

## 5.3. Mesons as non-thermal states

In the preceding section, we have understood the reason behind the lack of thermalization in the largeN c limit of (1+1)-d QCD. Our central result, that states containing heavy mesons do not thermalize in the long-time limit, directly implies a violation of the (strong) eigenstate thermalization hypothesis [18,19,56,57]: there must exist eigenstates within the spectrum where expectation values do not agree with the thermal average. This fits within a growing body of evidence [32-49] that system that exhibit confinement can exhibit an absence of thermalization, including some explicit examples that show violation of eigenstate thermalization [39,40].

In the language of the eigenstate thermalization hypothesis, the result from the previous sections can be reformulated as the fact that there exist non-thermal eigenstates, | ψ NT 〉 of the full Hamiltonian, H , of energy E , such that expectation values of local physical observables on this state do not agree with the thermal expectation value,

$$\langle \psi _ { \text{NT} } | O | \psi _ { \text{NT} } \rangle \neq \langle O \rangle _ { \beta },$$

chosen to have matching average energy 〈 H 〉 β = E . These states, | ψ NT 〉 , can be generated from the heavy-meson states of the previous section, as

$$| \psi _ { \text{NT} } \rangle = \lim _ { \tau \to \infty } e ^ { - \tau H } | n _ { 1 }, 0 ; n _ { 2 }, 0 ; \dots ; n _ { k }, 0 \rangle.$$

The process of taking the infinite-imaginary-time evolution of the heavy-meson states, amounts to selecting only the lowest lying eigenstate of H which has an overlap with the meson state. From the results of the previous section we can conclude that

$$| \psi _ { \text{NT} } \rangle = | n _ { 1 }, 0 ; n _ { 2 }, 0 ; \dots ; n _ { k }, 0 \rangle + \mathcal { O } \left ( \frac { 1 } { \sqrt { N _ { c } } } \right ),$$

that is, we have shown that eigenstates of the full Hamiltonian are well-approximated by the multi-meson states, being only weakly dressed by the interactions through O (1 / √ N c ) corrections. This is contrary to what is expected to happen in theories satisfying the eigenstate thermalization hypothesis, where any non-thermal eigenstates of the free Hamiltonian, H 0 should receive non-perturbative O (1) corrections once integrabilitybreaking interactions are included, rendering the eigenstates of the full Hamiltonian thermal.

The idea that theories with confinement and emergent mesonic excitations do not obey the eigenstate thermalization hypothesis is (in hindsight) not so surprising. Consider the limit case of a single, large meson, there is a priori no reason to expect that expectation values within such a state should agree with a thermal average at the same energy density, where many excitations exist in thermal equilibrium. An important aspect for such states to remain well-separated from the thermal continuum is that the decay of single, high-energy meson excitations should be suppressed. In the case considered here, this is achieved under a large N c prescription, whilst in the now well-studied case of the Ising model [36-40], this can easily be achieved with a weak longitudinal magnetic field.

There are numerous interesting directions to explore to further understand when, and why, mesons can be non-thermal. Tracing how the physics evolves away from the large N c limit would be interesting, although this requires non-perturbative methodologies. Hamiltonian truncation [97] may provide one such route, having already been applied to obtain low-energy eigenstates of largeN c QCD [98] and the NambuJonas-Lasino model, which also exhibits meson and baryon excitations [99]. One might expect that at intermediate values of N c and at short-to-intermediate time scales,

the physics still resembles the largeN c limit and statement made here apply to the prethermalization regime [83, 100-102]. This would be worthwhile studying. Another interesting regime to consider would be the deconfined phase, where it has been known that heavy mesons can persist [103]: Can these be interpreted in the framework of non-thermal states?

## 6. Conclusions

The largeN c limit of QCD has been well-studied since the seminal works of 't Hooft in the mid-1970s [66,70]. In light of recent works that have shown systems with confinement exhibit anomalous non-equilibrium dynamics and an absence of thermalization [32-49], we have return to this theory to address the question of whether the 't Hooft model thermalizes. Utilizing a powerful particle number/power of N c correspondence, and the annihilation pole structure of matrix elements within the theory, we have argued that the largeN c limit fails to thermalize. This statement holds even once non-perturbative corrections, which can arise at long times, are taken into account.

In particular, we have shown that simple initial states that contain few heavy meson excitations (with extensive total energy) do not thermalize. Here thermalize is used in the sense that expectation values of local observables do not relax, in the long time limit, to the relevant thermal result. At the heart of this physics is the local nature of meson excitations with respect to local observables, which implies there are no annihilation poles in matrix elements featuring one meson particle in both incoming and outgoing states. This absence of an annihilation pole prevents non-perturbative corrections that can bridge the O (1) gap between the t = 0 initial state expectation value and the thermal result. Thus, for this class of initial states, thermalization is avoided.

The absence of annihilation poles in one-particle-to-one-particle matrix elements is a crucial difference to analogous calculations in the Ising field theory and the sine-Gordon model [84-90]. There the corresponding matrix elements contain singularities that arise from the non-local soliton-like nature of excitations, and these lead to non-perturbative corrections in the time-evolution at long times. This allows for a O (1) difference between initial and final expectation values in these systems. Our work thus provides insight into a general mechanism in non-integrable (1+1)-d systems with confinement whereby thermalization can be avoided.

These results lend further credence to the burgeoning literature on the absence of thermalization in models with confinement [32-49]. We do note that our results and arguments apply to the largeN c limit of (1+1)-d QCD, and their applicability to the more interesting case of (3+1)-d and particle-collider phenomenology is not so evident. In the context of the quark-gluon plasma generated in heavy ion collisions, for example, it is expected that the collision energy is high enough to enter the deconfined phase. There is no obvious reason to expect the class of initial states that we have studied, containing few very heavy weakly-confined mesons, to have any relevance in that case. A more feasible future experimental testing ground could be quantum simulations of gauge

theories with cold atomic gases (see, e.g., refs. [104-106]), where engineering different non-equilibrium initial conditions and subsequent dynamics may be possible.

Additionally we have extended the concept of largeN c volume independence to nonequilibrium settings in confining gauge theories. These results are not limited to the lowdimensional 't Hooft model, but can be applied to more realistic gauge theories in higher dimensions. This opens a door for future studies of QCD non-equilibrium dynamics, by restricting oneself to evaluating only the leading corrections in the 1 /N c expansion, corresponding to highly non-diagonal matrix elements of physical observables.

To finish, we note that non-equilibrium dynamics can be strongly affected by the presence of few, special states (which may, for example, be thermodynamically unimportant). The confining Ising model is a nice example of this: the family of single meson eigenstates that sit above the thermal multimeson continuum are unimportant in equilibrium probes, such as the dynamical spin-spin correlation functions measured in inelastic neutron scattering. They are simply not observed in such probes, being overwhelmed by the proximate-in-energy thermal states (see, e.g., Ref. [60]). However, there is now plenty of evidence that such states can completely dominate non-equilibrium dynamics following simple quenches from simple, realisable initial states [36,39,40,48,49]. Thus when one finds interesting and unusual states, the question becomes: does my experiment realise and/or probe these? While current experiments on the quark-gluon plasma likely do not probe the states suggested here, perhaps their presence may be observed in different experimental set ups or systems in the future.

## Acknowledgments

We are grateful to Andrew James and Robert Konik for useful discussions and comments on the manuscript. A.C.C. is supported by the European Union's Horizon 2020 research and innovation programme under the Marie Skglyph[suppress] lodowska-Curie grant agreement 75009. N.J.R. is supported by the European Union's Horizon 2020 research and innovation programme under the Marie Skglyph[suppress] lodowska-Curie grant agreement 745944.

## References

- [1] Langen T, Geiger R and Schmiedmayer J 2015 Annu. Rev. Condens. Matter Phys. 6 201-217 URL https://doi.org/10.1146/annurev-conmatphys-031214-014548
- [2] Calabrese P, Essler F H L and Mussardo G 2016 J. Stat. Mech. 2016 064001 URL https: //doi.org/10.1088%2F1742-5468%2F2016%2F06%2F064001
- [3] Kinoshita T, Wenger T and Weiss D S 2006 Nature 440 900-903 URL http://dx.doi.org/10. 1038/nature04693
- [4] Polkovnikov A, Sengupta K, Silva A and Vengalattore M 2011 Rev. Mod. Phys. 83 (3) 863-883 URL https://link.aps.org/doi/10.1103/RevModPhys.83.863
- [5] Shuryak E 2005 Nucl. Phys. A 750 64 - 83 ISSN 0375-9474 quark-Gluon Plasma. New Discoveries at RHIC: Case for the Strongly Interacting Quark-Gluon Plasma. Contributions from the RBRC Workshop held May 14-15, 2004 URL http://www.sciencedirect.com/science/article/ pii/S0375947404011340

- [6] Blaizot J P, Gelis F, Liao J, McLerran L and Venugopalan R 2012 Nucl. Phys. A 873 68 -80 ISSN 0375-9474 URL http://www.sciencedirect.com/science/article/pii/ S0375947411006403
- [7] Pasechnik, Roman and Sumbera, Michal 2017 ˇ Universe 3 ISSN 2218-1997 URL http://www. mdpi.com/2218-1997/3/1/7
- [8] Shuryak E 2017 Rev. Mod. Phys. 89 (3) 035001 URL https://link.aps.org/doi/10.1103/ RevModPhys.89.035001
- [9] D'Alessio L, Kafri Y, Polkovnikov A and Rigol M 2016 Adv. Phys. 65 239-362 URL https: //doi.org/10.1080/00018732.2016.1198134
- [10] Rigol M, Dunjko V, Yurovsky V and Olshanii M 2007 Phys. Rev. Lett. 98 (5) 050405 URL https://link.aps.org/doi/10.1103/PhysRevLett.98.050405
- [11] Rigol M, Dunjko V and Olshanii M 2008 Nature (London) 452 854-858 ( Preprint 0708.1324 )
- [12] Korepin V E, Bogoliubov N M and Izergin A G 1993 Quantum Inverse Scattering Method and Correlation Functions Cambridge Monographs on Mathematical Physics (Cambridge University Press)
- [13] Nandkishore R and Huse D A 2015 Annu. Rev. Condens. Matter Phys. 6 15-38 URL https: //doi.org/10.1146/annurev-conmatphys-031214-014726
- [14] Alet F and Laflorencie N 2018 C. R. Phys. 19 498 -525 ISSN 1631-0705 quantum simulation / Simulation quantique URL http://www.sciencedirect.com/science/article/ pii/S163107051830032X
- [15] Parameswaran S A and Vasseur R 2018 Rep. Prog. Phys. 81 082501 URL https://doi.org/10. 1088%2F1361-6633%2Faac9ed
- [16] Vidmar L and Rigol M 2016 J. Stat. Mech. 2016 064007 URL https://doi.org/10.1088% 2F1742-5468%2F2016%2F06%2F064007
- [17] Ilievski E, Medenjak M, Prosen T and Zadnik L 2016 J. Stat. Mech. 2016 064008 URL https://doi.org/10.1088%2F1742-5468%2F2016%2F06%2F064008
- [18] Deutsch J M 1991 Phys. Rev. A 43 (4) 2046-2049 URL https://link.aps.org/doi/10.1103/ PhysRevA.43.2046
- [19] Srednicki M 1994 Phys. Rev. E 50 (2) 888-901 URL https://link.aps.org/doi/10.1103/ PhysRevE.50.888
- [20] Turner C J, Michailidis A A, Abanin D A, Serbyn M and Papi´ c Z 2018 Nat. Phys. 14 745-749 URL https://doi.org/10.1038/s41567-018-0137-5
- [21] Turner C J, Michailidis A A, Abanin D A, Serbyn M and Papi´ c Z 2018 Phys. Rev. B 98 (15) 155134 URL https://link.aps.org/doi/10.1103/PhysRevB.98.155134
- [22] Ho W W, Choi S, Pichler H and Lukin M D 2019 Phys. Rev. Lett. 122 (4) 040603 URL https://link.aps.org/doi/10.1103/PhysRevLett.122.040603
- [23] Khemani V, Laumann C R and Chandran A 2019 Phys. Rev. B 99 (16) 161101(R) URL https://link.aps.org/doi/10.1103/PhysRevB.99.161101
- [24] Lin C J and Motrunich O I 2019 Phys. Rev. Lett. 122 (17) 173401 URL https://link.aps.org/ doi/10.1103/PhysRevLett.122.173401
- [25] Choi S, Turner C J, Pichler H, Ho W W, Michailidis A A, Papi´ c Z, Serbyn M, Lukin M D and Abanin D A 2019 Phys. Rev. Lett. 122 (22) 220603 URL https://link.aps.org/doi/10. 1103/PhysRevLett.122.220603
- [26] Iadecola T, Schecter M and Xu S 2019 arXiv e-prints arXiv:1903.10517 ( Preprint 1903.10517 )
- [27] Bull K, Martin I and Papi´ c Z 2019 Phys. Rev. Lett. 123 (3) 030601 URL https://link.aps. org/doi/10.1103/PhysRevLett.123.030601
- [28] Moudgalya S, Rachel S, Bernevig B A and Regnault N 2018 Phys. Rev. B 98 (23) 235155 URL https://link.aps.org/doi/10.1103/PhysRevB.98.235155
- [29] Moudgalya S, Regnault N and Bernevig B A 2018 Phys. Rev. B 98 (23) 235156 URL https: //link.aps.org/doi/10.1103/PhysRevB.98.235156
- [30] Sala P, Rakovszky T, Verresen R, Knap M and Pollmann F 2019 arXiv e-prints arXiv:1904.04266

( Preprint 1904.04266 )

- [31] Khemani V and Nandkishore R 2019 arXiv e-prints arXiv:1904.04815 ( Preprint 1904.04815 )
- [32] Craps B, Kiritsis E, Rosen C, Taliotis A, Vanhoof J and Zhang H 2014 J. High Energy Phys. 2014 120 ISSN 1029-8479 URL https://doi.org/10.1007/JHEP02(2014)120
- [33] Craps B, Lindgren E and Taliotis A 2015 J. High Energy Phys. 2015 1-32 ISSN 1029-8479 URL https://doi.org/10.1007/JHEP12(2015)116
- [34] da Silva E, Lopez E, Mas J and Serantes A 2016 J. High Energy Phys. 2016 172 ISSN 1029-8479 URL https://doi.org/10.1007/JHEP06(2016)172
- [35] Myers R C, Rozali M and Way B 2017 J. Phys. A 50 494002 URL https://doi.org/10.1088% 2F1751-8121%2Faa927c
- [36] Kormos M, Collura M, Tak´ acs G and Calabrese P 2017 Nat. Phys. 13 246-249 URL http: //dx.doi.org/10.1038/nphys3934
- [37] Rakovszky T, Mesty´ an M, Collura M, Kormos M and Tak´ cs G 2016 Nucl. Phys. B 911 805 -845 ISSN 0550-3213 URL http://www.sciencedirect.com/science/article/pii/ S0550321316302541
- [38] H´ ods´gi K, Kormos M and Tak´cs G 2018 a a SciPost Phys. 5 (3) 27 URL https://scipost.org/ 10.21468/SciPostPhys.5.3.027
- [39] James A J A, Konik R M and Robinson N J 2019 Phys. Rev. Lett. 122 (13) 130603 URL https://link.aps.org/doi/10.1103/PhysRevLett.122.130603
- [40] Robinson N J, James A J A and Konik R M 2019 Phys. Rev. B 99 (19) 195108 URL https: //link.aps.org/doi/10.1103/PhysRevB.99.195108
- [41] Buyens B, Haegeman J, Van Acoleyen K, Verschelde H and Verstraete F 2014 Phys. Rev. Lett. 113 (9) 091601 URL https://link.aps.org/doi/10.1103/PhysRevLett.113.091601
- [42] Buyens B, Haegeman J, Verschelde H, Verstraete F and Van Acoleyen K 2016 Phys. Rev. X 6 (4) 041040 URL https://link.aps.org/doi/10.1103/PhysRevX.6.041040
- [43] Buyens B, Montangero S, Haegeman J, Verstraete F and Van Acoleyen K 2017 Phys. Rev. D 95 (9) 094509 URL https://link.aps.org/doi/10.1103/PhysRevD.95.094509
- [44] Buyens B, Haegeman J, Hebenstreit F, Verstraete F and Van Acoleyen K 2017 Phys. Rev. D 96 (11) 114501 URL https://link.aps.org/doi/10.1103/PhysRevD.96.114501
- [45] Ba˜ nuls M C, Cichy K, Cirac J I, Jansen K and K¨hn S 2017 u Phys. Rev. Lett. 118 (7) 071601 URL https://link.aps.org/doi/10.1103/PhysRevLett.118.071601
- [46] Akhtar A A, Nandkishore R M and Sondhi S L 2018 Phys. Rev. B 98 (11) 115109 URL https://link.aps.org/doi/10.1103/PhysRevB.98.115109
- [47] Park J, Kuno Y and Ichinose I 2019 arXiv e-prints arXiv:1903.07297 ( Preprint 1903.07297 )
- [48] Mazza P P, Perfetto G, Lerose A, Collura M and Gambassi A 2019 Phys. Rev. B 99 (18) 180302 URL https://link.aps.org/doi/10.1103/PhysRevB.99.180302
- [49] Lerose A, ˇ unkoviˇ Z c B, Silva A and Gambassi A 2019 Phys. Rev. B 99 (12) 121112 URL https://link.aps.org/doi/10.1103/PhysRevB.99.121112
- [50] Schecter M and Iadecola T 2019 arXiv e-prints arXiv:1906.10131 ( Preprint 1906.10131 )
- [51] Moudgalya S, Bernevig B A and Regnault N 2019 arXiv e-prints arXiv:1906.05292 ( Preprint 1906.05292 )
- [52] Bernien H, Schwartz S, Keesling A, Levine H, Omran A, Pichler H, Choi S, Zibrov A S, Endres M, Greiner M, Vuleti´ V and Lukin M D 2017 c Nature (London) 551 579-584
- [53] Essler F H L and Fagotti M 2016 J. Stat. Mech. 2016 064002 URL https://doi.org/10.1088% 2F1742-5468%2F2016%2F06%2F064002
- [54] Caux J S 2016 J. Stat. Mech. 2016 064006 URL https://doi.org/10.1088%2F1742-5468% 2F2016%2F06%2F064006
- [55] Lucini B and Panero M 2013 Phys. Rep. 526 93 -163 ISSN 0370-1573 URL http://www. sciencedirect.com/science/article/pii/S0370157313000021
- [56] Reimann P 2015 New J. Phys. 17 055025 URL https://doi.org/10.1088%2F1367-2630%2F17% 2F5%2F055025

- [57] Deutsch J M 2018 Rep. Prog. Phys. 81 082001 URL https://doi.org/10.1088%2F1361-6633% 2Faac9f1
- [58] Mondaini R and Rigol M 2017 Phys. Rev. E 96 (1) 012157 URL https://link.aps.org/doi/ 10.1103/PhysRevE.96.012157
- [59] McCoy B M and Wu T T 1978 Phys. Rev. D 18 (4) 1259-1267 URL https://link.aps.org/ doi/10.1103/PhysRevD.18.1259
- [60] Coldea R, Tennant D A, Wheeler E M, Wawrzynska E, Prabhakaran D, Telling M, Habicht K, Smeibidl P and Kiefer K 2010 Science 327 177-180 ISSN 0036-8075 URL http://science. sciencemag.org/content/327/5962/177
- [61] Morris C M, Vald´ es Aguilar R, Ghosh A, Koohpayeh S M, Krizan J, Cava R J, Tchernyshyov O, McQueen T M and Armitage N P 2014 Phys. Rev. Lett. 112 (13) 137403 URL http: //link.aps.org/doi/10.1103/PhysRevLett.112.137403
- [62] Wang Z, Schmidt M, Bera A K, Islam A T M N, Lake B, Loidl A and Deisenhofer J 2015 Phys. Rev. B 91 (14) 140404 URL http://link.aps.org/doi/10.1103/PhysRevB.91.140404
- [63] Wang Z, Wu J, Xu S, Yang W, Wu C, Bera A K, Islam A T M N, Lake B, Kamenskyi D, Gogoi P, Engelkamp H, Wang N, Deisenhofer J and Loidl A 2016 Phys. Rev. B 94 (12) 125130 URL https://link.aps.org/doi/10.1103/PhysRevB.94.125130
- [64] Bera A K, Lake B, Essler F H L, Vanderstraeten L, Hubig C, Schollw¨ck U, Islam A T M N, o Schneidewind A and Quintero-Castro D L 2017 Phys. Rev. B 96 (5) 054423 URL https: //link.aps.org/doi/10.1103/PhysRevB.96.054423
- [65] Fonseca P and Zamolodchikov A 2006 arXiv e-prints hep-th/0612304 ( Preprint hep-th/0612304 )
- [66] 't Hooft G 1974 Nucl. Phys. B 75 461 - 470 ISSN 0550-3213 URL http://www.sciencedirect. com/science/article/pii/0550321374900881
- [67] Eguchi T and Kawai H 1982 Phys. Rev. Lett. 48 (16) 1063-1066 URL https://link.aps.org/ doi/10.1103/PhysRevLett.48.1063
- [68] Kovtun P, ¨ nsal M and Yaffe L G 2007 U J. High Energy Phys. 2007 019 URL https://doi.org/ 10.1088%2F1126-6708%2F2007%2F06%2F019
- [69] Sulejmanpasic T, Shao H, Sandvik A W and Unsal M 2017 ¨ Phys. Rev. Lett. 119 (9) 091601 URL https://link.aps.org/doi/10.1103/PhysRevLett.119.091601
- [70] 't Hooft G 1974 Nucl. Phys. B 72 461 - 473 ISSN 0550-3213 URL http://www.sciencedirect. com/science/article/pii/0550321374901540
- [71] Callan C G, Coote N and Gross D J 1976 Phys. Rev. D 13 (6) 1649-1669 URL https: //link.aps.org/doi/10.1103/PhysRevD.13.1649
- [72] Einhorn M B 1976 Phys. Rev. D 14 (12) 3451-3471 URL https://link.aps.org/doi/10.1103/ PhysRevD.14.3451
- [73] Zhitnitsky A R 1985 Phys. Lett. B 165 405 -409 ISSN 0370-2693 URL http://www. sciencedirect.com/science/article/pii/0370269385912559
- [74] Li M 1986 Phys. Rev. D 34 (12) 3888-3893 URL https://link.aps.org/doi/10.1103/ PhysRevD.34.3888
- [75] Lenz F, Thies M, Levit S and Yazaki K 1991 Ann. Phys. (N.Y.) 208 1 - 89 ISSN 0003-4916 URL http://www.sciencedirect.com/science/article/pii/0003491691903426
- [76] Mermin N D and Wagner H 1966 Phys. Rev. Lett. 17 (22) 1133-1136 URL https://link.aps. org/doi/10.1103/PhysRevLett.17.1133
- [77] Hohenberg P C 1967 Phys. Rev. 158 (2) 383-386 URL https://link.aps.org/doi/10.1103/ PhysRev.158.383
- [78] Coleman S 1973 Commun. Math. Phys. 31 259-264 ISSN 1432-0916 URL https://doi.org/10. 1007/BF01646487
- [79] Burkardt M 1994 arXiv e-prints hep-ph/9409333 ( Preprint hep-ph/9409333 )
- [80] Manohar A V 1998 Large N QCD Probing the standard model of particle interactions. Proceedings, Summer School in Theoretical Physics, NATO Advanced Study Institute, 68th session, Les Houches, France, July 28-September 5, 1997. Pt. 1, 2 pp 1091-1169 ( Preprint hep-ph/

9802419 )

- [81] Hatta Y and Fukushima K 2004 Phys. Rev. D 69 097502 ( Preprint hep-ph/0307068 )
- [82] Cort´ es Cubero A 2016 J. Stat. Mech. 1608 083107 ( Preprint 1604.03879 )
- [83] Essler F H L, Kehrein S, Manmana S R and Robinson N J 2014 Phys. Rev. B 89 (16) 165104 URL https://link.aps.org/doi/10.1103/PhysRevB.89.165104
- [84] Schuricht D and Essler F H L 2012 J. Stat. Mech. 2012 P04017 URL https://doi.org/10. 1088%2F1742-5468%2F2012%2F04%2Fp04017
- [85] Calabrese P, Essler F H L and Fagotti M 2011 Phys. Rev. Lett. 106 (22) 227203 URL https: //link.aps.org/doi/10.1103/PhysRevLett.106.227203
- [86] Calabrese P, Essler F H L and Fagotti M 2012 J. Stat. Mech. 2012 P07016 URL https: //doi.org/10.1088%2F1742-5468%2F2012%2F07%2Fp07016
- [87] Calabrese P, Essler F H L and Fagotti M 2012 J. Stat. Mech. 2012 P07022 URL https: //doi.org/10.1088%2F1742-5468%2F2012%2F07%2Fp07022
- [88] Bertini B, Schuricht D and Essler F H L 2014 J. Stat. Mech. 2014 P10035 URL https: //doi.org/10.1088%2F1742-5468%2F2014%2F10%2Fp10035
- [89] Cort´ es Cubero A and Schuricht D 2017 J. Stat. Mech. 2017 103106 URL https://doi.org/10. 1088%2F1742-5468%2Faa8c2e
- [90] Horv´ ath D X, Kormos M and Tak´ acs G 2018 JHEP 2018 170 ISSN 1029-8479 URL https: //doi.org/10.1007/JHEP08(2018)170
- [91] Smirnov F A 1992 Form Factors in Completely Integrable Models of Quantum Field Theory (World Scientific) URL https://www.worldscientific.com/doi/abs/10.1142/1115
- [92] Mussardo G 2010 Statistical field theory (New York, NY: Oxford Univ. Press)
- [93] Delfino G, Mussardo G and Simonetti P 1996 Nuclear Physics B 473 469 - 508 ISSN 0550-3213 URL http://www.sciencedirect.com/science/article/pii/0550321396002659
- [94] Mussardo G 2010 Statistical Field Theory: An Introduction to Exactly Solved Models in Statistical Physics Oxford Graduate Texts (OUP Oxford) ISBN 9780199547586 see the discussion of Section 21.4, 'Non-locality and confinement'. URL https://books.google.nl/books?id= fakVDAAAQBAJ
- [95] Jaffe R and Mende P F 1992 Nucl. Phys. B 369 189 -218 ISSN 0550-3213 URL http: //www.sciencedirect.com/science/article/pii/055032139290383M
- [96] Babujian H, Fring A, Karowski M and Zapletal A 1999 Nucl. Phys. B 538 535 - 586 ISSN 05503213 URL http://www.sciencedirect.com/science/article/pii/S0550321398007378
- [97] James A J A, Konik R M, Lecheminant P, Robinson N J and Tsvelik A M 2018 Rep. Prog. Phys. 81 046002 URL https://doi.org/10.1088%2F1361-6633%2Faa91ea
- [98] Katz E, Marques Tavares G and Xu Y 2014 arXiv e-prints arXiv:1405.6727 ( Preprint 1405.6727 )
- [99] Azaria P, Konik R M, Lecheminant P, P´lmai T, Tak´cs G and Tsvelik A M 2016 a a Phys. Rev. D 94 (4) 045003 URL https://link.aps.org/doi/10.1103/PhysRevD.94.045003
- [100] Berges J, Bors´nyi a S and Wetterich C 2004 Phys. Rev. Lett. 93 (14) 142002 URL https: //link.aps.org/doi/10.1103/PhysRevLett.93.142002
- [101] Moeckel M and Kehrein S 2008 Phys. Rev. Lett. 100 (17) 175702 URL https://link.aps.org/ doi/10.1103/PhysRevLett.100.175702
- [102] Bertini B, Essler F H L, Groha S and Robinson N J 2015 Phys. Rev. Lett. 115 (18) 180601 URL https://link.aps.org/doi/10.1103/PhysRevLett.115.180601
- [103] Asakawa M and Hatsuda T 2004 Phys. Rev. Lett. 92 (1) 012001 URL https://link.aps.org/ doi/10.1103/PhysRevLett.92.012001
- [104] Zohar E, Cirac J I and Reznik B 2013 Phys. Rev. Lett. 110 (12) 125304 URL https://link.aps. org/doi/10.1103/PhysRevLett.110.125304
- [105] Surace F M, Mazza P P, Giudici G, Lerose A, Gambassi A and Dalmonte M 2019 arXiv e-prints arXiv:1902.09551 ( Preprint 1902.09551 )
- [106] Celi A, Vermersch B, Viyuela O, Pichler H, Lukin M D and Zoller P 2019 arXiv e-prints arXiv:1907.03311 ( Preprint 1907.03311 )
Axel Cortés Cubero Neil J. Robinson
Lack of thermalization in (1+1)-d QCD at large $N_c$
Motivated by recent works aimed at understanding the status of equilibration
and the eigenstate thermalization hypothesis in theories with confinement, we
return to the 't Hooft model, the large-$N_c$ limit of (1+1)-d quantum
chromodynamics. This limit has been studied extensively since its inception in
the mid-1970s, with various exact results being known, such as the quark and
meson propagators, the quark-antiquark interaction vertex, and the meson decay
amplitude. We then argue this model is an ideal laboratory to study
non-equilibrium phenomena, since it is manifestly non-integrable, yet one
retains a high level of analytic control through large-$N_c$ diagrammatics. We
first elucidate what are the non-equilibrium manifestations of the phenomenon
of large-$N_c$ volume independence. We then find that within the confined
phase, there is a class of initial states that lead to a violation of the
eigenstate thermalization hypothesis, i.e. the system never thermalizes. This
is due to the existence of heavy mesons with an extensive amount of energy, a
phenomenon that has been numerically observed recently in the quantum Ising
chain.