Thermalization of Yang-Mills theory in a $(3+1)$ dimensional small lattice system

We study the real-time evolution of SU($2$) Yang-Mills theory in a $(3+1)$ dimensional small lattice system after interaction quench. We numerically solve the Schr{\"o}dinger equation with the Kogut-Susskind Hamiltonian in the physical Hilbert space obtained by solving Gauss law constraints. We observe the thermalization of a Wilson loop to the canonical state; the relaxation time is insensitive to the coupling strength, and estimated as $\tau_{\rm eq}\sim 2\pi/T$ with temperatures $T$ at steady states. We also compute the vacuum persistence probability (the Loschmidt echo) to understand the relaxation from the dynamics of the wave function.


I. INTRODUCTION
How an isolated quantum system reaches thermal equilibrium is one of the fundamental problems in modern physics.In particular, the thermalization of a nonabelian gauge theory is important for understanding the nature of quark-gluon plasmas observed in relativistic heavyion collision experiments (See Refs.[1,2] for recent reviews).The analysis of relativistic heavy-ion collision experiments using a hydrodynamic model implies that hydrodynamics can be applied from a very early stage after the collision (∼ 0.5 fm/c).However, a microscopic calculation based on kinetic theory shows that the time scale of thermalization with a small QCD coupling α S is of order α −13/5 S Q −1 s , where Q s is the characteristic momentum scale of gluons inside the colliding nuclei.This time scale is orders of magnitude larger than the one expected in hydrodynamic models [3,4].On the other hand, analyses based on the gauge/gravity duality imply that the thermalization time is the order of the inverse of temperature 1/(πT ) in the large colors and large 't Hooft coupling limit [5][6][7][8].This rapid thermalization has been thought of as the universal property of the strongly coupled gauge theories, and the quark-gluon plasmas produced in relativistic heavy-ion collision experiments are thought of as strongly coupled, while the conventional plasmas are weakly coupled.
Hydrodynamic behavior has also been observed in small systems of pp and pA collisions [9][10][11].This observation is unexpected because it is usually thought that thermalization does occur due to a large number of degrees of freedom.There are two possibilities: One is that this is the property of quantum field theories.Since the dimension of the Hilbert space is infinite, thermalization does occur even in small systems.The other is that hydrodynamics works well before thermalization where the pressure is isotropic, which is called hydrodynamization.This possibility has been intensively studied in recent years (see Refs. [2,12] for recent reviews).
In this paper, we discuss the first possibility of quantum thermalization.For this purpose, we study the thermalization of the SU(2) Yang-Mills theory in a small isolated system.We employ the Kogut-Suskind Hamiltonian formulation on a single cubic lattice with open boundary conditions [13] and numerically solve the Schrödinger equation.To mimic the situation in heavyion collision experiments, we employ interaction quench from the strong gauge coupling limit to a weak gauge coupling.This enables us to avoid the aforementioned problem of whether thermalization, hydrodynamization, or isotropization.Furthermore, the thermalization mechanism will be purely quantum because there is no kinetic regime on the small lattice.
The advantage of the Hamiltonian formulation is free from the so-called sign problem in Monte Carlo simulations, that is, the difficulty of importance sampling due to the complexity of the path-integral weight [14].On the other hand, the disadvantage is the exponentially large Hilbert space.However, as shown below, we overcome this difficulty by considering a small system and explicitly solving the Gauss law constraints, which numerously reduces the size of physical Hilbert space, and enables us to access the real-time dynamics of the Yang-Mills theory using the standard classical computers.Our finding has a substantial impact on developing fields of classical or quantum simulations of lattice gauge theories (see Refs. [15,16] for review).This paper is organized as follows.In Sec.II, we briefly review the Kogut-Susskind Hamiltonian formulation and show how to construct the physical space and the matrix element of Hamiltonian.In Sec.III, we discuss the realtime dynamics by solving the Schrödinger equation.We show that the results of our numerical simulations exhibit thermalization, and the thermalization time is of order of the Boltzmann time, τ eq ∼ 2π/T .Section IV is devoted to summary and outlook.In Appendix A, we show the eigenvalues of the Hamiltonian near the ground state.(2) Yang-Mills theory.The gauge field Uµ(x) is defined on a link.The chromo-electric fields EL and ER are defined on the ends of a link, and satisfy the Gauss law constraint on the vertices, e.g., ER6 + ER9 + EL10 = 0 at the vertex indicated by the black dot.The minimal Wilson loop operator is defined as the product of the link operators on the edges of a unit plaquette colored in the figure.

II. HAMILTONIAN FORMULATION A. formulation
We review the Hamiltonian formulation of the lattice SU(2) Yang-Mills theory, which is often referred to as the Kogut-Susskind Hamiltonian formulation [13].We consider a cubic lattice (See Fig. 1).The gauge fields U µ (x) are defined on a link emanating from a site x and terminating at x+ê µ , with êµ=x,y,z being the unit vector along the µ direction.U µ (x) is a 2 × 2 matrix-valued operator, and we can apply the local SU(2) transformation to U µ (x) from left-/right-hand side.Using the generators E a L (x, µ), and E a R (x, µ) [a = x, y, z], which are nothing but the left-or right-chromoelectric fields defined on the ends of a link, the SU(2) algebras are represented as where σ a=x,y,z are the Pauli matrices, and other commutation relations vanish.abc is the Levi-Civita symbol with xyz = 1.In the Kogut-Suskind Hamiltonian formulation, the generators are not independent but related through the constraint =: E 2 (x, µ). ( The Hamiltonian is given as the sum of electric and magnetic parts, H = H E + H B , with + (h.c.), (7) where P is the set of plaquettes, and "h.c." represents the Hermitian conjugate.K is the coupling constant, which is inversely proportional to the square of the gauge coupling g.Therefore, we use the words, strong and weak coupling, for small and large K, respectively.H E is the electric part of the Hamiltonian, which has the same form as continuum theory.H B is the lattice version of the magnetic part of the Hamiltonian; it involves a nonlocal but gauge-invariant operator of U µ (x), which is the famous Wilson loop operator (see Fig. 1).The Schrödinger equation with the Hamiltonian defines the dynamics on the physical Hilbert space |Ψ that satisfies the Gauss law constraints: The Gauss law constraints simply state that the total electric field at a site x must vanish.
For numerical implementation, we rewrite the Kogut-Susskind Hamiltonian using the so-called Schwinger bosons [17].In SU (2), the Hamiltonian of the chromoelectric fields (6) is the same as that of the quantum rotor.Therefore, the electric field operator can be understood as the angular momentum operator, and represented by using the creation and annihilation operators of the spin doublet bosons (Schwinger bosons) as where a i=↑,↓ , a † i=↑,↓ (b i=↑,↓ , and b † i=↑,↓ ) are the annihilation and creation operators of the Schwinger bosons, which are defined on the left (right) end of a link.In terms of the Schwinger bosons, the constraint (5) implies where i b i are the number operators of the Schwinger bosons.Therefore, the total number of Schwinger bosons living on the edges of a link must be the same.Using the Schwinger boson representation, the electric part of the Hamiltonian is written only by N L or N R .
Next, we consider the magnetic part of the Hamiltonian.It is known that the link operator can be written using the Schwinger bosons as Using the commutation relations between creation and annihilation operators , Eqs. ( 9), ( 10), ( 12) and ( 13) reproduce those between E L,R , and U in Eqs. ( 1)-( 4).We can label the Hilbert space of the gauge theory by the number of eigenvectors of the harmonic oscillators.This enables us to understand the complex wave function of the gauge theory from the simple picture of the occupation dynamics of bosons.Furthermore, by truncating the max occupation number of the Schwinger bosons, we can obtain the finite-dimensional Hilbert space with manifestly keeping the gauge symmetry.
Two remarks are in order: (I) Truncating the Schwinger boson occupation numbers at a certain value is equivalent to representing the Kogut-Susskind theory with the j-dimensional irreducible representation of SU(2) up to j max = (N L + N R )/2.The limit j max → ∞ may recover Wilson's formulation of lattice gauge theories.(II) The magnetic part of the Hamiltonian (7) does change the number of the Schwinger bosons on each link with satisfying the constraint (5), while the electric part of the Hamiltonian (6) just counts their numbers.Therefore, they can be understood as the kinetic and interaction terms of the Schwinger bosons.Without magnetic interactions (K = 0) corresponding to the strong coupling limit, the gauge theory is reduced to free harmonic oscillators.As K increases, the fluctuations by the magnetic Hamiltonian become relevant, and then the gauge theory becomes strongly correlated.In what follows, we study such strongly correlated dynamics of the Yang-Mills theory by quenching the magnetic Hamiltonian and solving the time-dependent Schrödinger equation after the quench.We note that larger K demands larger j max for full quantitative analysis since the effect of the truncation of the Hilbert space becomes more relevant as the typical occupation number increases.
The Gauss law constraint (8) needs to be satisfied at each vertex.
Since j i has no upper bound, the dimension of the physical Hilbert space is infinite.In numerical simulations, we truncate the spin's maximum value, j max , to make the dimension of the Hilbert space finite.We show the j max dependence of the dimension of the physical Hilbert space in Table I.

C. Matrix element of Hamiltonian
Let us evaluate the matrix element of the Hamiltonian H = H E + H M .Because |j is an eigenstate of H E , the matrix element of H E is just the sum of eigenvalues of spins: For the magnetic part, more calculations are involved.Since H M consists of the sum of plaquettes, let us, first, focus on the single plaquette tr(U Here, we employed the cyclic property of the trace.It is useful to express the matrices on each vertex as where we define This expression enables us to express the plaquette as the sum of L's, Each L s s (b, a) locally acts on the Hilbert space, so that it is enough to consider the action of L s s (b, a) on a single vertex.Let |j 1 , j 2 , j 6 be a local state on the vertex (1, 2, 6) ∈ V .Since L s s (b, a) is written by the creation and annihilation operators, it is easy to calculate the action of L s s (b, a).
Using these definitions, we find the matrix element of the magnetic part as where s p is a vector whose components are defined as In summary, the physical state is given in Eq. ( 16), and the matrix element of the Hamiltonian is expressed as × (sgn(i, j))

III. REAL-TIME SIMULATION A. Interaction quench
We consider the single cubic lattice with open boundary conditions, shown in Fig. 1.We truncate the Schwinger boson occupation number at For example, if we consider the lowest truncation j max = 1/2, the dimension of the local Hilbert space is 5.
In the number basis |N L↑ N L↓ |N R↑ N R↓ , these are explicitly given as |00 |00 , |10 |10 , |10 |01 , |01 |10 , |01 |01 .Therefore, the dimension of the full Hilbert space is 5 12 ∼ 0.2 billion.The full Hilbert space is so large that we cannot manage it in numerical simulations except the lowest j max .However, the majority of the Hilbert space represents the redundancy associated with the gauge symmetry, and we need only the subspace (physical Hilbert   space) obtained by solving the Gauss law constraints (8).This process significantly reduces the dimension of the Hilbert space, e.g., from 5 12 to 32 for j max = 1/2.By explicitly solving the Gauss law constraints, we can do numerical simulations with larger j max .
We numerically solve the time-dependent Schrödinger equation, i∂ t |Ψ(t) = H|Ψ(t) , in the physics Hilbert space and Hamiltonian constructed in the previous section.As an initial state, we choose the Fock vacuum defined by a i (x, µ)|Ψ(0) = b i (x, µ)|Ψ(0) = 0.The Fock vacuum is the eigenstate of the electric Hamiltonian (6), that is, the ground state of the Hamiltonian at the strong coupling limit K = 0. We study the real-time dynamics after the magnetic Hamiltonian is switched on at t = 0. We solved the time-dependent Schrödinger equation based on the leap-frog type discretization.We decompose the time-dependent Schrödinger equation into two real-valued equations: We regard the real and imaginary parts as "position" and "velocity" and apply the leap-frog integrator.This method is applicable only when the Hamiltonian is realvalued in some basis 1 .The numerical resources needed to obtain the following results, e.g., with j max = 4 (d = 87, 426, 119) are 262 TFlops*hr for each K.

B. Thermalization time
We show the time evolution of the Wilson loop after the interaction quench in Fig. 3.We clearly see the Wilson loop rapidly reaches some equilibration value and fluctuates around it.The fluctuation around the  II.
long time average decreases as K increases, which is less than 1% for K > 5, while it is about 15% for K = 1 (see Table II).To reveal whether the Wilson loop reaches the thermal state or not, we computed the canonical ensemble average of the same Wilson loop operator, (Tr e −βH trU )/Z, where Tr represents the trace over the physical Hilbert space, β = 1/T is the inverse temperature, and Z = Tre −βH is the partition function.We choose β so that the canonical average of the Hamiltonian (Tre −βH H)/Z equals the expectation value in real-time evolution Ψ(t)|H|Ψ(t) = Ψ(0)|H|Ψ(0) = E, where E is the total energy that is equal to zero in our state.We found that the long-time average is in accord with the canonical ensemble average.This implies the Wilson loop gets equilibrated to the thermal state.Since the fluctuations for K = 1, 2 are not small, we focus on the time evolutions for K = 5, 10, 15, 25 in the following.
Let us evaluate the time scale of relaxation to the thermal equilibrium.We show the log plot of the deviation from the long-time average in Fig. 4. The log plot shows linear decreasing behavior in time, which implies the exponential damping of the deviation, trU (t) − trU ∼ e −t/τeq , with trU being the time average after thermalization.Although the strength of the expectation values strongly depends on K, the time scale of thermalization is insensitive.There is an ambiguity to determine the thermalization time due to oscillations of expectation < l a t e x i t s h a 1 _ b a s e 6 4 = " < l a t e x i t s h a 1 _ b a s e 6 4 = " E h S O G r y D p 0 f E V j r X B 2 s / F h w k r h c = " > A A A E J H i c h Z P P a x N R E M e / 6 f q j 1 h 9 t 9 S J 4 C c a K i I S r c w W u 8 w a 6 1 Z 3 2 x v l r 7 w 9 C J U s 5 c w T / L + v Y H y P / n P A = = < / l a t e x i t > K = 5 < l a t e x i t s h a 1 _ b a s e 6 4 = " a c j q s c l D C 6 Q M q E p 6 6 K + i E E 2 h i / Y = " > A A A E J X i c h Z P P a x N R E M e / 6 f q j 1 h / 9 4 U X w E o w V E Q l v p d A i C A U v Q i 9 t 2 p h C r W V 3 8 9 I 8 u r 9 4 + x J b l x y 9 e P D a g i c F B f H P 8 O L B o y 3 0 5 F k 8 V v D i w d n J l q q h m 7 d s d u a 7 8 5 l 5 M 9 n n x r 5 K j B C H p R H r z N l z 5 0 c v j F 2 8 d P n K + M T k 1 O M k 6 m h P 1 r 3 I j / S q 6 y T S x f E n O 2 v z r F P 2 i 5 k m x b T o 6 n P H / 5 H m i s 2 h t O K v Y J A c X l f z t / v 3 P G q s 9 O j 0 2 P + f l U G j c a 9 q z 1 R t e 2 m m M n 8 z P 0 i j u I 4 b u E 2 H Z R b z e I R F 1 K l w G 6 + w i z 3 r v f X F + m o d 9 E N H S j l z F f 8 s 6 9 s f l Z H n c g = = < / l a t e x i t > K = 10 < l a t e x i t s h a 1 _ b a s e 6 4 = " y F z L 2 y t j B B S S / l v 4 H / 8 / V P z O N J A = " > A A A E J X i c h Z P P a x N R E M e / 6 f q j r T / a 6 q X g J R g r I h L e l k q L I B S 8 C F 7 a t D G F W s v u 5 q V 5 d H / x 9 i V a l x y 9 e P B a w Z O C g v h n e P H g U Y W e P I v H F n r R 9 h 9 8 B 7 t y x r r z l 7 S M w T 7 j n g + i F N O S W 9 S x n b l L m X e 9 l M j 5 Q i z p z I P c b N I e w q U 5 m l S S + z X 0 w 8 z e t k s 0 7 Z K 4 4 / 3 l m b Z 5 2 y X 8 w 0 K a Z F V 5 8 7 + o 8 0 V 2 w O p R V / B Y P k 8 L q a v 9 2 / 5 1 F j p U e n x / 7 / r A w a j d m q P V e 1 7 e W 5 y u K 1 / C C N 4 g q u 4 g Y d l n k s 4 j 6 W U K f C b b z E L l 5 Z 7 6 0 v 1 l f r e z 9 0 p J Q z l / H P s n 7 8 A a h H 5 3 c = < / l a t e x i t > K = 15 < l a t e x i t s h a 1 _ b a s e 6 4 = " + R 9 5 B X E U 6 9 0 8 E + D B Q l T / 7 + 9 p m g I = " > A A A E J X i c h Z P P a x N R E M e / 6 f q j r T / a     II.
values.We, here, employ a linear fit using the peaks of the normalized deviation log | trU (t) / trU − 1| from t = 0 to t = 25β in Fig. 4. The time scales of thermalization are summarized in Table II, which are typically τ eq ∼ 6.5 × β ∼ 2πβ = 2π/T .The time scale 2π/T is known as the Boltzmann time [23,24].
For typical temperature of the quark-gluon-plasma produced in RHIC, T = 200 MeV, τ eq ∼ 2π/T is 6 fm/c.This is an order of magnitude larger than the time scale expected from the hydrodynamic model (∼ 0.5 fm/c), and the thermalization time scale 1/(πT ) ∼ 0.3fm/c, observed in calculations based on the gauge/gravity duality [5][6][7][8].
Next, to deepen understanding of the thermalization from the dynamics of the wave function, we compute the vacuum persistent probability, which is also known as the Loschmidt echo or fidelity in the context of quantum chaos [25] and dynamical quantum phase transition [26].The Loschmidt echo is defined as

and quantifies the deviation of the state < l a t e x i t s h a 1 _ b a s e 6 4 = " q
r u H g r J P M w y 7 6 z j 6 z H v r I j 9 o O d X p k r o h z x X v b w a c M j v L 1 + B h 7 U J t / e X f 8 9 l P X w q a F 1 Q a X u X E M D l m n H A j s I S I l 7 c f p 8 5 / X 7 3 v r j 0 m x 0 n x 2 y n 9 j F A T t h X 7 A P 2 f n l f F j j p X 3 K L p F 5 S T 1 7 V F / i l C l y k w 4 y z 6 y M 6 z i A z t l x 1 i H 7 P 2 0 P 6 3 y y n v y L p G z R z W 7 F F 9 i l 0 P E e + i x g 5 7 7 s R b 1 9 A J J 4 u l r e T s w P 4 K 7 T q l y k w 4 y z 6 y M 6 z i A z t l x 1 i H 7 P 2 0 P 6 3 y y n v y L p G z R z W 7 F F 9 i l 0 P E e + i x g 5 7 7 s R b 1 9 A J J 4 u l r e T s w P 4 K 7 T q e 7 u N V X S D 6 i X P Z x t e A 5 P m 7 f A / f r E + 8 f r f 0 Z y n V x 1 d C + Z i V m r q E J i 5 S x w A p 8 Q q J a 7 D 6 / e 3 B 0 v v a y P B P O s k / s F 1 b x k Z 2 x E 6 x D d n / b n 1 d 5 + Q N 5 l 8 h 5 R z W 7 F F 9 i l 0 P E u + i x j Z 5 7 s R b 1 9 A p J 4 u l b e d s w N 4 S 7 T q x I U o j n S U 9 m 7 M < l a t e x i t s h a 1 _ b a s e 6 4 = " Z S W r N 6 P h h G K z u 5 2 2 I 7 u z y + y 0 g p s e v X A n H j x h 4 o H w E z y q i X f j g Z 9 g P G J i T D z 4 9 n U J a M N 2 J r P z 3 j f v e 2 / e 2 3 m W 7 4 h A M 3 a W G k l n b t y 8 N T q W v T 1 + 5 2 4 u P z H 5 L P A 6 y u Y V 2 3 M 8 t W W Z A X e E 5 B U t t M O 3 f M V N 1 3 J 4 1 d p 9 G p 1 X u 1 w F w p O b + s D n O 6 7 Z k q I p b F M j V M / P 1 y z e E j L k e x 1 C 5 n r Z F / W w p t y C a + 7 3 n i x m a 1 w 2 r p z W 8 0 x A e j J j P 4 4 T 1 T o k L d n + 8 m Z t q n V I e j K n g T Z N n H 3 e x T 9 S F L E x l C 3 o F Q w y h 8 d V 9 H a v 1 q N M S N Q 9 x v + 9 M i h U F 0 r G Y s k w 1 h e K K w / j R h q F K X g A s 9 g s S 7 A C q 7 A G F Q x 8 B B / g E 3 x O f 0 3 / z q Q y 6 b 7 p S C r m 3 I N / R m b 8 L / K i 9 1 Q = < / l a t e x i t > j max = 4 FIG.6.Time evolution of the Wilson loop defined on the bottom plaquette in Fig. 1 in units of β for K = 10, and jmax = 3/2 (yellow), jmax = 2 (cyan), jmax = 5/2 (magenta), jmax = 3 (red), jmax = 7/2 (green), and jmax = 4 (blue).The values of β for K = 10, and jmax = 3/2 (yellow), jmax = 2 (cyan), jmax = 5/2 (magenta), jmax = 3 (red), jmax = 7/2 (green), and jmax = 4 (blue) are 0.03034, 0.04093, 0.05127, 0.06173, 0.07134, and 0.07996, respectively.The dashed lines show the canonical ensemble averages.at time t [|Ψ(t) ] from the initial state t = 0 [|Ψ(0) ].We show the time evolution of the logarithm of the Loschmidt echo after the interaction quench in Fig. 5.The wave function rapidly spreads out to the entire physical Hilbert space, and then after the time scale where the Wilson loop gets equilibrated, the spreading also stops and fluctuates around equilibrium values.We can see three characteristic time regions: early, intermediate, and late time.At the early time, the logarithm of the Loschmidt echo shows a quadratically decreasing, log P vac (t) −η 2 t 2 .The value of η depends on the strength of the interaction.We can nicely fit the data as η ∼ 2.5 × K, which is independent of the temperature.At the intermediate time, log P vac (t) is linearly damping with oscillations, −t/τ LE .Again, this is oscillating, so that we employ a linear fit using points at the peak positions.It is, however, not easy to evaluate τ LE for K ≤ 5 because no clear peaks are found.We here only evaluate τ LE for K > 5, which leads to the typical time scale, τ LE ∼ 0.3 × β ∼ 1/(πT ).This is comparable to the thermalization time scale observed in calculations based on the gauge/gravity duality [5][6][7][8].At the late time, log P vac (t) fluctuates around −10.

C. jmax dependence
Finally to see the dependence of the truncation of spins (i.e., the j max dependence), we show the time evolution of the Wilson loop by changing j max in Fig. 6 with an intermediate coupling K = 10.We also show the j max dependence of the normalized deviation of the same Wilson loop from its long-time average in Fig. 7.We see that the relaxation time scale is insensitive to the choice l y k w 4 y z 6 y M 6 z i A z t l x 1 i H 7 P 2 0 P 6 3 y y n v y L p G z R z W 7 F F 9 i l 0 P E e + i x g 5 7 7 s R b 1 9 A J J 4 u l r e T s w P 4 K 7 T q l y k w 4 y z 6 y M 6 z i A z t l x 1 i H 7 P 2 0 P 6 3 y y n v y L p G z R z W 7 F F 9 i l 0 P E e + i x g 5 7 7 s R b 1 9 A J J 4 u l r e T s w P 4 K 7 T q < l a t e x i t s h a 1 _ b a s e 6 4 = " Z S W r N 6 P h 9 L A J v E z Y Q 5 b e S W B V 2 W 4 = " > A A A E U n i c h V N B T x N B F H 6 l V a G i L X g x 4 d J Y I Y S Y M k t I I C Y m J F 4 4 Q q G W h G K z u 5 2 2 I 7 u z y + y 0 g p s e v X A n H j x h 4 o H w E z y q i X f j g Z 9 g P G J i T D z 4 9 n U J a M N 2 J r P z 3 j f v e 2 / e 2 3 m W 7 4 h A M 3 a W G k l n b t y 8 N T q W v T 1 + 5 2 4 u P z H 5 L P A 6 y u Y V 2 3 M 8 t W W Z A X e E 5 B U t t M O 3 f M V N 1 3 J 4 1 d p 9 G p 1 X u 1 w F w p O b + s D n O 6 7 Z k q I p b F M j V M / P 1 y z e E j L k e x 1 C 5 n r Z F / W w p t y C a + 7 3 n i x m a 1 w 2 r p z W 8 0 V W Y j Q K g 4 I R C 0 W I x 5 o 3 k X o N N W i A B z Z 0 w A U O E j T K D p g Q 4 N w G A x j 4 i O 1 A i J h C S d A 5 h x 5 k k d t B K 4 4 W J q K 7 + G 2 h t h 2 j E v X I Z 0 B s G 6 M 4 u B Q y C z D N v r E T d s 6 + s F P 2 n f 2 5 1 l d I P q K 7 H O B u w S N c b t 8 D 9 + u 5 w / s b v 4 Z y X d w 1 t C 9 Z i T f X 0 I R l u r H A D H x C o l z s P r / 7 6 s 3 5 x u P y d D j D 3 r E f m M U x O 2 M f M Q / Z / W m / X + f l t + R d I u c l 5 e x S f I l V D h H v o s c 2 e u 7 F W l T T C y S J p 6 / l P Y e 5 I d x N Y k W S Q r x A e j J j P 4 4 T 1 T o k L d n + 8 m Z t q n V I e j K n g T Z N n H 3 e x T 9 S F L E x l C 3 o F Q w y h 8 d V 9 H a v 1 q N M S N Q 9 x v + 9 M i h U F 0 r G Y s k w 1 h e K K w / j R h q F K X g A s 9 g s S 7 A C q 7 A G F Q x 8 B B / g E 3 x O f 0 3 / z q Q y 6 b 7 p S C r m 3 I N / R m b 8 L / K i 9 1 Q = < / l a t e x i t > of j max , in particular for j max > 2, where the dimension of the Hilbert space is larger than 10 6 .In this case, the quantitative discussion of the relaxation time is possible within the numerically reachable j max .
To see the j max dependence of the absolute value of the Wilson loop.We show the canonical average of the Wilson loop (corresponding to stationary values in Fig. 6) by changing j max in Fig. 8 with an intermediate coupling K = 10.Using the linear and quadratic polynomial fitting, we estimate the canonical average with j max → ∞ as 1.33 ± 0.08 (linear), and 1.50 ± 0.13 (quadratic).The largest three j max are used for the linear fitting, while all j max are used in the quadratic fitting.Here the errors are estimated from the 95% confidence interval.We show the fitting curves in Fig. 8.Although the two estimations give the consistent results within the error bars, each data strongly depends on j max as seen in Fig. 8, and results in large extrapolation errors.This result implies that we may need larger j max for the complete quantitative research.

IV. SUMMARY AND OUTLOOK
We have studied the real-time evolution of the SU(2) Yang-Mills theory in a (3 + 1) dimensional small lattice system after interaction quench.We have numerically solved the Schrödinger equation in the reduced Hilbert space obtained by explicitly solving the Gauss law constraints.We have observed the thermalization to the canonical state; the relaxation time τ eq is insensitive to the strength of the coupling constant, and scaled by the Boltzmann time 2π/T .The observed thermalization is very rapid compared with conventional matters, but it is still an order of magnitude larger than the one expected from the hydrodynamic model.We hope that our numerical simulations in small systems share essential features of nonequilibrium dynamics with real QCD, although we need to confirm it by conducting a more comprehensive study in future works, e.g., checking the j max , system size, and initial-state dependences of the relaxation time, changing the lattice geometry, generalizing to the SU(3) group, and so on.In particular, the strong j max dependence is observed in the absolute value of the Wilson loop, although the relaxation time is less insensitive to j max .We may elaborate on these in future research.
Furthermore, using our formulation, we can attack important problems of nonequilibrium QCD.For example, we can compute the Kubo formula and estimate transport coefficients in a small system.We can also compute the so-called out of time-order correlators, and confirm whether the lattice Yang-Mills theory saturates the maximum bound on the quantum Lyapunov exponent conjectured on the basis of the gauge/gravity duality [27].

FIG. 1 .
FIG. 1. Schematic representation of the lattice SU(2) Yang-Mills theory.The gauge field Uµ(x) is defined on a link.The chromo-electric fields EL and ER are defined on the ends of a link, and satisfy the Gauss law constraint on the vertices, e.g., ER6 + ER9 + EL10 = 0 at the vertex indicated by the black dot.The minimal Wilson loop operator is defined as the product of the link operators on the edges of a unit plaquette colored in the figure.

12 FIG. 2 .
FIG.2.Labels of links on the single cubic lattice model.
H / P D l 2 4 P N + 7 V F d Y N 8 I D 9 x F + / J P v m C + + D D X / 7 H D V p 7 p 7 N z Z J 7 r P Y e 6 P s c u K 9 S H m L G P m d P c y 3 o 6 V o o 4 e S z 3 F F a m s H 3 p d I i C A U v g p c 2 a U y h 1 r K 7 e U k e 3 V + 8 f Y n W J U c v g l d 7 8 K T Q g / h n e B E 8 S g + 9 e B e P F b x 4 c H a y p W r o 5 i 2 b n f n u f G b e T P a 5 s a 8 S I 8 R B a c I 6 d f r M 2 c l z U + c v X L w 0 6 y o g 5 i V r B e v z 3 e f 7 x 4 t 3 6 9 N p 7 f E W / G D u n g j D s U n 6 i P s / v T e L c n a a 8 4 e E v O M e w 6 4 f k h T T k n v U s Y 2 Z e 7 l X j b T Y 6 W I M 6 d y T 3 F n C L v C V G Z p 0 s v s F x P b e Z 1 s 1 i l 7
L H D w p f r T s 0 A V S 4 F X g 7 y 7 A d O + h d 7 R N w 5 P X w v L + U x r l l 2 y Z / J / w R 7 Z H d 3 A b r z p V 5 s 8 f 4 4 Y f Y D 8 8 7 m 7 g b K Q l X N Z W d 7 M p V d W w 7 + I Y g o z m K M H X 8 Q K 1 r G B Y m D v B G c 4 j 7 x I S W l S m u 6 k S p F Q M 4 5 v I W U + A b t B j J A = < / l a t e x i t > t/ < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 8 y 5 s w 9 d 8 7 M 6 I 4 p P Z / o O a a 0 t L a 1 d 3 R 2 x b t 7 e v s S y f 6 B L c 8 u u 4 b I G b Z p u 9 u 6 5 g l T
e 7 f O 6 6 8 p u 0 T m O e 0 5 o P o S u 5 y g P s C M P c y c 5 l 7 W 0 2 O l i N N n c s 9 g a Q K 7 S V R m K d R N 8 o u J F 3 m d r N c J e c X x J y v r U a 8 T 8 o u Z N s Z 0 8 B p x x 9 9 I U c X 2 R F r Q X z B O T q 6 r 6 N 8 9 3 Y 8 6 K d n p s f 4 / K + N G c 7 l m r d Q s a 3 2 5 u v o w P 0 g z c A N u w R 0 8 L P d g F Z 7 A G j S w 8 A f 4 A k P 4 N v 3 d m D W u G w u j 0 K l S z l y D f 4 Z x 8 y + / L A A R < / l a t e x i t > htrU ⇤ i < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 Z e g 2 f p u o r c e o R D 3 y G R D b x i g O L o X M P E y z n 2 y P n b B D t s 9 + s b/ X + g r J R 3 S X X d w t m M P l 9 j x w v z b 5 / u 7 K 6 U C u i 7 u G 1 g U r 8 e Y a G v C U b i w w A 5 + Q K B e 7 x + + 8 / X i y 8 q w 0 H T 5 k X 9 h v z O I z O 2 b f M Q / Z + W N / X e a l T + R d I u c N 5 e x S f I l V D h H v o M c W e u 7G W l T T c y S J p 6 / l b c L s A O 4 q s S J J I Z 4 n P Z m x E 8 e J a h 2 S l m x / c b M W 1 T o k P Z l T R 5 s G z h 7 v / B 8 p i l g f y B b 0 C v q Z g + M q e r u X 6 1 E i J O o e 4 2 q v 9 A u V h a L x u G g Y y w u F x Q d x I 4 3 C P b g P M 9 g s T 2 A R X s I S l D H w B / g G P + A g f Z Q + y 6 Q y 6 Z 7 p U C r m 3 I H / R m b 8 H + 7 U 9 1 M = < / l a t e x i t > j max = 3 < l a t e x i t s h a 1 _ b a s e 6 4 = " f e G t n f y 5 2 2 A 7 u z 6 + y 0 A p s e u f g F O H j S x I P x M 3 A i U b 8 A B z 6 C 8 Y i J B z n w 9 n U J a M N 2 J r v z 3 m / e 7 / 3 b m b V 8 R w S a s b P U S D p z 5 + 7 o 2 L 3 s / Q c P x y d y k 1 r E T d s 6 + s F P 2 n f 2 5 1 l d I P q K 7 H O B u w S N c b t 8 D 9 + u 5 w / s b v 4 Z y X d w 1 t C 9 Z i T f X 0 I R l u r H A D H x C o l z s P r / 7 6 s 3 5 x u P y d D j D 3 r E f m M U x O 2 M f M Q / Z / W m / X + f l t + R d I u c l 5 e x S f I l V D h H v o s c 2 e u 7 F W l T T C y S J p 6 / l P Y e 5 I d x N Y k W S Q r

<
l a t e x i t s h a 1 _ b a s e 6 4 = "V 3 E S k 9 L d R T n 0 T x r 2 4 M 1 B i 4 j r v m o = " > A A A C a 3 i c h V G 7 S g N B F D 1 Z X z E + E m M j a h E M E a s 4 K w H F S r S x 1 M S 4 A Q 1 h d 5 3 o k n 2 x O w n E 4 A / Y W Q l a K Y i I n 2 H j D 1 j 4 C a K d g o 2 F N 5 s F 0 a D e Y W b O n L n n z p k Z z T U N X z D 2 G J F 6 e v v6 B 6 K D s a H h k d F 4 Y i y 5 7 T t 1 T + d F 3 T E d r 6 S p P j c N m x e F I U x e c j 2 u W p r J F a 2 2 1 t 5 X G t z z D c f e E k 2 X l y 1 1 3 z

4 F
X g 7 y 7 A d O + h d 7 R N w 5 P X w v L + U x r l l 2 y Z / J / w R 7 Z H d 3 A b r z p V 5 s 8 f 4 4 Y f Y D 8 8 7 m 7 g b K Q l X N Z W d 7 M p V d W w 7 + I Y g o z m K M H X 8 Q K 1 r G B Y m D v B G c 4 j 7 x I S W l S m u 6 k S p F Q M 4 5 v I W U + A b t B j J A = < / l a t e x i t > t/< l a t e x i t s h a 1 _ b a s e 6 4 = " H / P D l 2 4 P N + 7 V F d Y N 8 I D 9 x F + / J P v m C + + D D X / 7 H D V p 7 p 7 N z Z J 7 r P Y e 6 P s c u K 9 S H m L G P m d P c y 3 o 6 V o o 4 e S z 3 F F a m s H e g 2 f p u o r c e o R D 3 y G R D b x i g O L o X M P E y z n 2 y P n b B D t s 9 + s b / X + g r J R 3 S X X d w t m M P l 9 j x w v z b 5 / u 7 K 6 U C u i 7 u G 1 g U r 8 e Y a G v C U b i w w A 5 + Q K B e 7 x + + 8 / X i y 8 q w 0 HT 5 k X 9 h v z O I z O 2 b f M Q / Z + W N / X e a l T + R d I u c N 5 e x S f I l V D h H v o M c W e u 7 G W l T T c y S J p 6 / l b c L s A O 4 q s S J J I Z 4 n P Z m x E 8 e J a h 2 S l m x / c b M W 1 T o k P Z l T R 5 s G z h 7 v / B 8 p i l g f y B b 0 C v q Z g + M q e r u X 6 1 E i J O o e 4 2 q v 9 A u V h a L x u G g Y y w u F x Q d x I 4 3 C P b g P M 9 g s T 2 A R X s I S l D H w B / g G P + A g f Z Q + y 6 Q y 6 Z 7 p U C r m 3 I H / R m b 8 H + 7 U 9 1 M = < / l a t e x i t > j max = 3 < l a t e x i t s h a 1 _ b a s e 6 4 = " f e G t n f y i o T f i m i 9 S E t f c 3 p l G t W Y = " > A A A E V H i c h V R P T x N B F H + l V b A q F L i Y e G m s E E N M m W 1 M I C Y k J F 4 8 Q q G W h E K z u 5 2 2 A 7 u z 6 + y 0 A p s e u f g F O H j S x I P x M 3 A i U b 8 A B z 6 C 8 Y i J B z n w 9 n U J a M N 2 J r v z 3 m / e 7 / 3 b m b V 8 R w S a s b P U S D p z 5 + 7 o 2 L 3 s / Q c P x y d y k 1 N v A q + j b F 6 x P c d T G 5 Y Z c E d I X t F C O 3 z D V 9 x 0 L Y d X r d 1 X 0 X 6 1 y 1 U g P L m u 9 3 2 + 5 Z o t K Z r C N j V C 9 Z x R s 3 h L y J C / 7 R A y 1 8 v u 1 M O a c v O u u d d b W p g v Z W t c N m 7 s 1 3 M F V m Q 0 8 o O C E Q s F i M e K N 5 k 6 h B o 0 w A M b O u A C B w k a Z Q d M C H B u g g E M f M S 2 I E R M o S R o n 0 M P s s j t o B V H C x P R X X y 3 U N u M U Y l 6 5 D M g t o 1 R H H w U M v M w w 0 7 Z F 3 b O f r C v 7 C e 7 u N V X S D 6 i X P Z x t e A 5 P m 7 f A / f r E + 8 f r f 0 Z y n V x 1 d C + Z i V m r q E J i 5 S x w A p 8 Q q J a7 D 6 / e 3 B 0 v v a y P B P O s k / s F 1 b x k Z 2 x E 6 x D d n / b n 1 d 5 + Q N 5 l 8 h 5 R z W 7 F F 9 i l 0 P E u + i x j Z 5 7 s R b 1 9 A p J 4 u l b e d s w N 4 S 7 T qx I U o j n S U 9 m 7 M V x o l 6 H p C X b X 2 f W p l 6 H p C d z G m j T x N n n X X 0 j R R E b Q 9 m C T s E g c 3 h c R W f 3 Z j / K h E S 3x / j / r g w K 1 V L R e F E 0 j N V S Y f l p f J H G 4 D E 8 g W d 4 W R Z g G V 7 D C l Q w 8 B E c w z f 4 n j 5 N / 8 3 g X 6 J v O p K K O d P w z 8 i M X w L W L f f M < / l a t e x i t > j max = 7/2

TABLE I .
jmax dependence of the dimension of the physical Hilbert space d.

TABLE II .
Time average after thermalization trU , the canonical ensemble average trU can of the Wilson loop, and the corresponding inverse temperature β. σ, τeq, and τLE are the normalized standard deviation ( trU / trU − 1) 2 , the relaxation time to the thermal state, and the time scale of the exponential decay of the Loschmidt echo, respectively.