Thermalization of long range Ising model in different dynamical regimes: a full counting statistics approach

We study thermalization of transverse field Ising chain with power law decaying interaction $\sim 1/r^{\alpha}$ following a global quantum quench of the transverse field to two different dynamical regimes. We quantify the thermalization behavior by comparing the full probability distribution function (PDF) of the evolving states with the corresponding thermal state given by the Gibbs canonical ensemble (GCE). To this end, we use matrix product state (MPS) based time dependent variational principle (TDVP) algorithm to simulate both real time evolution following a global quantum quench and the finite temperature density operator. We observe that thermalization is strongly suppressed in the region with strong confinement for all the interaction strength $\alpha$ considered whereas thermalization occurs in the region with weak confinement.


Introduction
The investigation of non-equilibrium dynamics in isolated many-body systems has garnered significant attention in recent decades, owing to advancements in the manipulation of synthetic quantum systems in laboratory settings [1][2][3][4][5][6][7][8][9][10][11] and the development of analytical and numerical techniques [12][13][14][15][16][17][18][19][20][21][22][23][24][25].An enduring question in quantum many-body dynamics pertains to the potential thermalization of a closed system that has been perturbed from equilibrium.Thermalization implies that the long-term behavior of a dynamical system can be anticipated using the principles of statistical mechanics.Generally, in the case of a non-integrable closed system, one would expect thermalization in accordance with the Eigenstate Thermalization Hypothesis (ETH) [26][27][28][29][30]. Nonetheless, certain studies have presented contradictory evidence, at least within their specific regime and time scales of investigation [31][32][33][34].In a scenario where a closed system is initially prepared in a generic state, denoted as |ψ i 〉 (which is not an eigenstate of the Hamiltonian), and subsequently evolved unitarily with a non-integrable Hamiltonian, Ĥ, thermalization is said to occur if the local observables eventually relax to an equilibrium state that corresponds to the predictions of the thermal ensemble, 〈 Ô〉 t|t→∞ → 〈 Ô〉 eq = 〈 Ô〉 CGE = Tr( ρβ Ô) Equation (1) indicates thermalization according to the microcanonical ensemble (MCE).The summation encompasses all the eigenstates of the Hamiltonian within a narrow energy range δE centered around the initial energy E i = 〈ψ i | Ĥ|ψ i 〉.The normalization factor N E i ,δE tallies the energy eigenstates within this range, spanning 2δE.This approach requires either full diagonalization [28,35,36] or partial diagonalization centered on the initial energy density [37] of the Hamiltonian.Consequently, computational limitations arise as a function of the system size.Equation (2) implies thermalization in accordance with the canonical Gibbs ensemble (CGE).The trace is performed over the density operator ρβ , defined as the inverse temperature β = 1/T , which is determined by the system's initial energy E i = Tr( ρβ Ĥ) Tr( ρβ ) .Further elaboration on how to extract β and ρβ is provided in section 3.2.Equations ( 1) and (2) are recognized as the conditions for strong thermalization.An alternative weak thermalization condition occurs when the time-averaged local order parameter converges to the thermal prediction [33,34].
Disordered systems demonstrating many-body localization impede thermalization [38][39][40][41].In clean systems, dynamical confinement hinders information propagation and the thermalization process [42][43][44][45][46][47][48][49][50][51][52][53].The Long Range Ising model (LRIM) exhibits confinement [44,45] as outlined in Appendix D. A recent study [4] observed the suppression of thermalization in the confined regime of LRIM simulated with trapped ions.Furthermore, employing high-scale exact diagonalization, the validity of ETH has been examined for various interaction strength parameters, denoted as α (see section 2), in LRIM.The results indicate that a strong ETH typically holds at least within the range α ≥ 0.6 [35].In our previous work [54], we explored the relaxation of order parameter statistics following a global quench, revealing two distinct dynamical regimes based on the gaussification of the full counting statistics (FCS) of subsystem magnetization.Building on this foundation, the present study investigates the thermalization of LRIM under the CGE framework following a global quench into different dynamical regimes.We evaluate thermalization using the most rigorous criteria by comparing the FCS of the time-evolving state post-global quench with that of the corresponding thermal state.
The remainder of this paper is organized as follows: In Section 2, we introduce the model, the order parameter, its distribution, and a metric for quantifying the proximity of the timeevolved state to the thermal state.Section 3 provides comprehensive details of the numerical methods, quench protocol, and extraction of effective temperatures associated with the global quench.Specific details on simulating finite temperature density operators and calculating the full counting statistics of the order parameter is provided in Appendix B.1 and Appendix B.2 respectively.Section 4 presents the outcomes of our study.Appendix B.3 details the error analysis of the numerical results and Appendices C and D provide supplementary information on thermal phase transitions and correlation propagation in the Long Range Ising model.Finally, we conclude by summarizing our findings and suggesting potential avenues for future research in Section 5.

Model and Methods
We investigate the ferromagnetic long-range Ising model (LRIM) described by the Hamiltonian in Equation (3), where ŝµ i , µ = x, y, z is the spin one-half matrices at site i.We consider open boundary condition, that is relevant to existing experimental setups.The ferromagnetic interaction between two spins falls as the inverse power of the distance between them and is parameterized by the interaction strength parameter α.For α ≤ 1, the inverse power-law interaction series diverges with the lattice size and is normalized using the Kac normalization constant, as defined in Equation (4).
This normalization ensures the intensivity of the energy density in the regime α ≤ 1.The static and dynamic behaviors of this model are strongly influenced by the interaction strength parameter α.At α = ∞, the model simplifies to the transverse field Ising model (TFIM), which can be solved exactly by mapping it to a system of spinless fermions through Jordan-Wigner transformations [55].TFIM exhibits a quantum phase transition from the ferromagnetic phase to the paramagnetic phase at h = |J|/2.This quantum phase transition persists as α decreases, with the transition point shifting towards higher values of the magnetic field h [5,56,57].At the opposite extreme of α = 0, we have a fully connected regime that is amenable to analytical treatment for both the static and dynamic properties [54,58,59].For α < 2, this model displays a long-range ferromagnetic order at low finite temperatures [60,61].Given the absence of spontaneous symmetry breaking in finite systems and the 2 symmetry of the model in Equation (3), the finite-temperature states with ferromagnetic order also exhibit 2 symmetry (see appendix C).The regime with α < 2 is particularly intriguing and features a wealth of exotic phenomena such as prethermalization [11], nonlinear propagation of light cones [62,63], dynamical phase transitions [2,3,54,[64][65][66][67][68][69], and dynamical confinement [4,44,45].Furthermore, this model has garnered significant attention owing to its experimental relevance, particularly in systems involving trapped ions with adjustable transverse field strengths and interaction ranges [1][2][3][4]11].
The complete information of a generic time evolving quantum state, expanded in the computational basis |ψ t 〉 = {σ i } C {σ i } (t) |σ 1 , σ 2 , . . ., σ i , . . ., σ N 〉, is encapsulated within the set of time dependent coefficients {C {σ i } (t)}.In many-body systems, these coefficients scales exponentially with the number of spins, rendering their study exceedingly challenging.A common approach for investigating dynamics in such systems is to monitor the evolution of the expectation value of a local observable 〈ψ t | Ô|ψ t 〉, such as the order parameter in systems exhibiting order-disorder transitions.A more robust strategy involves tracking the full probability distribution function (PDF) of this observable, which provides comprehensive information on quantum fluctuations in the system, including all moments and cumulants.Specifically, when the operator Ô is diagonal in computational basis, the corresponding PDF is defined as which represents the histogram of the squared coefficients of the many-body wave function within the range of possible outcomes of the measurements of Ô.The PDF allows for straightforward calculation of moments (or cumulants) of any order.In this study, we employ the eigenvectors of the total spin operator in the longitudinal direction, that is, Ŝ x = N i=1 ŝx i , as our computational basis.In this context, the order parameter of interest is the longitudinal magnetization defined for a subsystem of size l within a system of N spins: The observable Ml is suitable because it typically relaxes to a stationary state [70], ultimately approaching a stationary statistical distribution in a subsystem of dimension l.In addition, because Ml is diagonal in computational basis, the definition in Equation ( 5) applies.The probability distribution function of the subsystem magnetization Ml within a generic state ρt (whether pure or mixed) is given by which can be Fourier-transformed into an integral form: where G l (θ , t) = Tr ρt e iθ Ml denotes the moment-generating function.Given that the Hamiltonian (3) involves a system of spin one-half particles, the values of m span either integers or half-integers within the range m ∈ − l 2 , − l 2 + 1, . . ., l 2 − 1, l 2 depending on whether l is even or odd.In Appendix B.2 we illustrate the detailed calculation of G l (θ , t) with matrix product state (MPS) representation.Historically, the PDF has been studied as the full counting statistic (FCS) of electron fluctuations in mesoscopic systems [71][72][73].More recently, FCS has been explored in quantum many-body systems in both equilibrium and non-equilibrium scenarios [18,19,43,54,[74][75][76][77][78][79].
To assess thermalization, we introduce a metric called "Distance to Thermalization", DT(t), initially introduced in [80].This metric quantifies the Euclidean distance between the probability distribution function (PDF) of the order parameter at time t following a quantum quench, denoted as P l (m, t), and the corresponding thermal PDF, represented as P TH l (m).Mathematically, it is defined as It is noteworthy that the convergence of the PDF provides a more rigorous criterion for thermalization than the convergence of the expectation value.This is because the former implies the latter, whereas the reverse is not necessarily true.A similar approach has been employed in previous studies to investigate thermalization dynamics [37,[79][80][81].Comprehensive details of how to extract the thermal state corresponding to a global quantum quench are discussed in Section 3.2.In cases where the system undergoes thermalization, DT(t) is expected to converge to zero in the long time limit.

Real and imaginary time evolution
The numerical simulations in this study are classified into two distinct categories: • real time evolution of pure state following a global quench.
• simulation of finite temperature density operators.
For both of these simulation tasks, we employ the MPS-based Time Dependent Variational Principle (TDVP) algorithm [23,24] with second order integration scheme.This choice affords us a significant advantage over exact diagonalization methods, allowing us to simulate systems of much larger sizes than can be accommodated by the current exact diagonalization techniques.
The quench protocol implemented in this study, with energy rescaled to |J| = 1, is as follows: At time t = 0, the system is prepared in the ground state of the Hamiltonian Ĥi (α, 0), which is a 2 symmetric Greenberger-Horne-Zeilinger (GHZ) state oriented along the longitudinal direction: GHZ state is characterized by the PDF which is sharply bimodal with two peaks at m = l/2 and m = −l/2 respectively.Equation ( 10) can be explicitly represented as an exact Matrix Product State (MPS) with bond dimension χ = 2. Subsequently, a global quench is initiated along the transverse field h to a final Hamiltonian Ĥ f (α, h f ) and the system is evolved unitarily using the expression |ψ t+d t 〉 = e −id t Ĥ f |ψ t 〉.The evolution is monitored by calculating the Full Counting Statistics (FCS) of the order parameter at each time step.The details of the quench protocol is pictorially represented in Figure 1 (a).The finite temperature density operator is simulated by an imaginary time evolution starting from a maximally mixed state at infinite temperature.Additional details pertaining to the calculation of the thermal density operator are provided in Appendix B.1.For both sets of simulations, we maintain a fixed maximum bond dimension of the MPS at χ max = 128.Furthermore, Trotter time steps of d t = 0.05 and dβ = 0.001 are used for real and imaginary time evolution, respectively.There is a finite time-step error of O(d t 3 ) per time step and O(d t 2 ) per unit time [82].In Appendix B.3 we access the accuracy of the TDVP data by comparing the TDVP results with the exact results obtained by the full diagonalization of a system of size N = 14.Furthermore, we test the convergence of the data by calculating the relative error for three increasing bond dimensions.

Extraction of effective temperature of a global quench
A global quantum quench Ĥi (α, 0) − → Ĥ f (α, h) in an isolated system adds an extensive amount of energy to the system.Consequently, the system relaxes to a state at a higher energy level than the ground state of the post-quench Hamiltonian [70], where |ψ 0 〉 is the ground state of the post-quench Hamiltonian Ĥ f (α, h).The left hand side of Equation ( 11) is a conserved quantity because the real time evolution of |ψ t 〉 is unitary.
For every global quantum quench we can attribute an effective temperature β eff which is the temperature at which the thermal energy density above the ground state of the post-quench Hamiltonian matches the conserved energy density of the system, The effective temperature is extracted by solving Equation (12).The left hand side of the equation is trivially calculated as and the right hand side can be calculated for a series of β by numerically solving Equation ( 23) and calculating the energy density at each instance.The precision of β eff depends on the trotter steps dβ in the solution to equation (23).In Fig. 2 we plot the numerical solution of equation (12).The energy density attributed to quench (represented by the black dashed line) in our setup is independent of the post-quench parameters because the spin-spin interaction term in the Hamiltonian (3) is normalized with the Kac normalization (4), whereas the expectation value h 〈ψ i | j ŝz j |ψ i 〉 taken over the transverse field term is trivially zero.If we extend the simulation to a larger β (i.e., lower temperature), all curves will converge to the ground state energy density of Ĥ f at the corresponding post-quench parameters.Once β eff is extracted, we can calculate the corresponding thermal PDF, P TH (m) = P β eff (m), using equation (8).≈ 0.5 for α ≤ 2 [64,85].When h f < h dyn c , the system is in the dynamical ferromagnetic phase, which strongly retains the ferromagnetic order of the initial GHZ state following a global quantum quench.This is evident from the persistent oscillation of P l (m, t) around P GHZ l (m).Conversely, when h f > h dyn c , the system transits to the dynamical paramagnetic phase, characterized by the rapid dissolution of the initial ferromagnetic order of the initial GHZ state following a global quantum quench.This is signified by the Gaussification of P l (m, t) [54,64].The comprehensive dynamical phase diagram and universality behavior related to the dynamical phase transition in LRIM remain active areas of investigation [66,[86][87][88].Our primary objective is to examine the convergence behavior of the metric DT(t) in two distinct dynamical phases of the LRIM.The convergence of DT(t) towards zero is an indicator of thermalization within a particular phase under consideration.We initialize the system as crucial as this region exhibits an interesting landscape encompassing both finite temperature phase transitions [89,90] and dynamical phase transitions [2,3,54,[64][65][66][67][68][69].Specifically, we consider three distinct values of interaction strength, namely α = 0.0, 1.5, 1.9.At α = 0.0, the system exhibits integrability because of its full connectivity and complete permutation symmetry, thereby leading to a lack of thermalization [17].On the other hand, the choices of α = 1.5 and α = 1.9 are motivated by the relatively faster equilibration and Gaussification of the PDF following a quench in the dynamical paramagnetic phase, as previously observed [54].

Quench to dynamical ferromagnetic regime
Figure 3 shows the temporal evolution of DT(t) following a global quantum quench of the transverse field to h f = 0.3 with α = 0.0, 1.5, 1.9 for subsystem sizes l = 20, 60, 100.Notably, all these points belong to the dynamical ferromagnetic phase [54,64].For all three quenches, a persistent oscillation in DT(t) is evident, indicating that the initial ferromagnetic order is strongly retained and thermalization is suppressed.This behavior aligns with the relaxation mode represented by Path 3 in figure 1(b).Specifically, α = 0.0 is in the integrable regime; therefore, thermalization is expected to be absent [26] whereas we anticipate thermalization for quenches with α = 1.5 and α = 1.9.The apparent suppression of thermalization can be attributed to the confinement behavior.The long-range interaction of the model effectively confines low-energy domain wall kinks into heavier quasiparticles that typically travel slower than free quasiparticles, thereby suppressing the spread of correlations in the system [44,45].Consequently, thermalization is still expected but only at significantly longer time scales [47].Appendix D details the confinement behavior in LRIM where we observe the spreading of connected correlation function ŝx k ŝx k+∆ c for α = 1.9 and h f = 0.3 shows a strong temporal suppression.In figure 3(d), (e), and (f), the colored scattered plots depict P l (m) as a function Figure 4: First row: Time evolution of the metric DT(t) following a global quantum quench to three interaction strength values α ∈ {0.0, 1.5, 1.9} and transverse field h f = 0.6 at three different subsystem sizes l = {20, 60, 100}.All three points are in dynamical paramagnetic phases [54,64].Second row: P l (m) versus m for m ∈ [0, l/2] with l = 100 at four different time slices t = {2, 6, 20, 50}.P l (m) versus m for m ∈ [−l/2, 0) is its mirror image.The black dashed curve represents the thermal PDF, P TH (m) attributed to the corresponding global quantum quenches. of m at four distinct time slices.The black dashed curve represents the Probability Density Function (PDF) of the expected thermal state, P TH l (m).We observe that the time-evolving P l (m) oscillates persistently around P TH l (m).Of particular importance is the observation that, in all three cases, the thermal PDFs are bimodal, indicating the presence of long-range ferromagnetic order.This observation suggests that if the system eventually thermalizes for these post-quench parameters at extended time scales, it would exhibit a long-range ferromagnetic order.This finding further strengthens the argument that this is indeed a dynamical ferromagnetic phase.

Quench to dynamical paramagnetic regime
Figure 4 illustrates the temporal evolution of DT(t) following a global quantum quench of the transverse field to h f = 0.6, with α = 0.0, 1.5, 1.9 for subsystem sizes l = 20, 60, 100.These points are located within the dynamical paramagnetic phase [54,64].Notably, these quenches exhibit a distinct relaxation behavior of DT(t) compared with the previous cases.In Figure 4(a), we observe rapid equilibration for all values of l.However, it is essential to highlight that DT(t) remains at or above the order of O(10 −1 ) following equilibration, which suggests a lack of thermalization.This behavior aligns with expectations, because α = 0 represents an integrable point.For α = 1.5, DT(t) does not exhibit stable equilibration (see figure 4 (b)); Finally, when α = 1.9, we observe equilibration for l = 60, 100 (see figure 4 (c)).DT(t) exhibits a stable oscillation around a constant value of approximately O(10 −3 ).A more comprehensive picture is shown in Fig. 4(f), where the late-time PDF perfectly overlaps with the corresponding thermal PDF represented by a black dashed curve.This is indicative of thermalization of the corresponding quench.Although we observed signatures of thermalization, the system is not in a de-confined phase [4,91].In Appendix D, the connected correlation function ŝx for α = 1.9 and h f = 0.6 still exhibits weaker temporal suppression.A recent study observed a de-confinement transition for a system of up to 31 spins for a much higher value of the transverse field [4].This suggests that while strong confinement suppresses thermalization, signatures of thermalization can still be detected in the presence of weak confinement.To further support this observation we study the post quench temporal evolution of domain wall kinks defined as, k counts the number of nearest neighbor kinks in the x direction within subsystem l.Because confinement bounds the domain walls kinks into heavier quasiparticles, it is a relevant parameter to study.In Figure 5, we illustrate the temporal evolution of the average domain-wall kinks, 〈 k〉, following a global quantum quench.As anticipated, quenches to the dynamical ferromagnetic phase with h f = 0.3 display persistent oscillations around the thermal value, indicating a lack of thermalization.Conversely, for quenches to the dynamical paramagnetic phase with h f = 0.6, we observe distinctly different post-quench behavior.In the case of α = 0, the domain wall kinks equilibrate to a stable value that differs from the expected thermal value, as expected because it is an integrable point.This observation complements the post-quench behavior of DT, as depicted in Figure 4 (a).Although α = 1.5 is a non-integrable point, thermalization is not observed within the simulation time.At later times, a stable prethermal plateau, close but distinct from the expected thermal value, becomes apparent.Conversely, for a quench corresponding to α = 1.9, the average domain wall kinks converge to the expected thermal value.Notably, before reaching the thermal value, the kink density exhibits a relatively stable prethermal plateau until time t ≃ 35.This relaxation mode, which is characterized by two time scales, is represented by Path 2 in Figure 1(b).This discovery provides another robust indicator of thermalization in weakly confined regimes.

Conclusion
We investigate the relaxation dynamics of the long-range Ising model subsequent to a global quantum quench of the transverse field, assessing the thermalization on a computationally viable time scale according to the canonical Gibbs ensemble (CGE).The model is non-integrable for all values of α except at the extremes (α = 0.0, ∞), where we anticipate thermalization following a global quantum quench.However, the long-range Ising model exhibits confinement, which suppresses correlation spreading and eventually suppresses thermalization.Starting from the Greenberger-Horne-Zeilinger (GHZ) state, we quench the system to two distinct dynamical regimes.As anticipated, robust confinement suppresses thermalization for smaller quenches, specifically in the dynamical ferromagnetic region, where the metric DT exhibits persistent oscillations characteristic of the masses of bound mesons.Conversely, for quenches to the dynamical paramagnetic region, a notably different behavior emerges.The persistent oscillation diminishes and the DT relaxes more rapidly.Although conclusive evidence of thermalization for α = 1.5 is not observed within the simulation time, compelling indications of the thermalization surface for α = 1.9 are based on the relaxation of DT.This observation gains additional support from the convergence of the domain wall kinks to the expected thermal value.
The time evolved state is given by |ψ t 〉 = e −i Ĥ t |ψ 0 〉, where Ĥ is the post-quench Hamiltonian 3. We proceed by expanding |ψ 0 〉 in the eigenbasis, of the post-quench Hamiltonian, where q j = E j ψ 0 .Further expanding E j in the computational basis, E j = n c j n |n〉, we can derive the expression for q j as, The post-quench state is where X n (t) = N j=0 q j c j n e −i E j t .We can now calculate the time evolution of the expectation value of a generic parameter Ô as, If |n〉 is the simultaneous eigenket of the order parameter Ô then 18 becomes, Finally, with the full eigenvalues of hamiltonian at hand we can also calculate the energy density corresponding to a thermal density matrix ρβ ,

B Simulations details
In this section we present the details of numerical simulation complementary to the results in the main text.

B.1 Simulation of finite temperature density operator
The finite temperature states can be simulated by casting the density operator as locally purified tensors [92,93].The thermal density operator is defined by Gibbs distribution where β = 1 T is the inverse temperature.At β = 0 (infinite temperature) the state is maximally mixed and is given as the tensor product of local identities ρ0 and d is the dimension of the local Hilbert space (for spin 1 2 , d = 2).The density operator for any finite temperature (non-zero β) is We keep the density operator operator in locally purified form ρ = † at each stage where is represented as tensor where σ i = d, k i = d are the physical index and the Kraus index, they remain fixed, and 1 ≤ c i ≤ χ ma x is the bond index and χ max is the maximum value of bond dimension.The density operator initialized at infinite temperature can now be purified to a finite temperature in trotterized steps Equation ( 23) can be simulated using imaginary time TDVP (−id t → −dβ) in only the half section of the density operator operator and never contracting the X and X † layer during the evolution, thus strictly preserving the locally purified form.Figure 6 shows the tensor notation of the infinite temperature density operator ρβ=0 which is a tensor product of identity matrices of size (d, d), where d is the physical dimension.Rather than working with the density operator as an MPO we represent the density operator in the locally purified form [93,94] which is positive semi-definite by construction and keep it in locally purified form at every stage of the thermal purification process.In Fig. 7 we represent ρβ=0 in the locally purified form β=0 † β=0 , where the index in purple is an auxiliary index called the Krauss index.
we can now evolve one of the halves ( or † ) as shown in equation ( 23) and the evolution on the other half is its trivial conjugate.This approach is computationally efficient as we can work with cheaper MPS instead of more expensive MPDO.In Fig. 8 one half of the ρβ=0 in locally purified form is shown, form here on we will only work with this half.Algebraically, β=0 can be written as For the system of spin one-half particles we choose X as as shown in Fig. 9.This particular choice is taken to preserve the trace of the density operator, Figure 9: Choice of X σ i ,k i to preserve the trace of ρ.
Finally, we reshape β=0 from a string of 2 × 2 matrices to a string of four legged tensors of shape (1, 2, 2, 1) as shown in Fig. 10, which is a MPS of bond dimension 1.Now that we have our initial state as an MPS, we can simulate a finite temperature density operator by solving the equation ( 23), Numerically, equation ( 23) can be solved for long-range spin systems through imaginary time evolution, (where id t is transformed into dβ) using the Time-Dependent Variational Principle (TDVP).The TDVP algorithm employed for simulating the thermal state remains fundamentally identical to that used for the real-time evolution of the pure state, with the distinction of an additional auxiliary Krauss index.However, in the thermal purification of a closed system, the Krauss index becomes obsolete because all physical operators act solely on the physical index, and the Krauss indices are contracted among themselves [93].Figure 11 illustrates the tensor network diagram for computing the expectation value of a two point operator Ôi Ôj acting on site i and j within the thermal state ρβ .

B.2 Calculating full counting statistics with MPS
The Central object in the calculating the full probability distribution function of an order parameter is the moment generating function G l (θ , t) = Tr( ρt e iθ Ml ).For pure state, the density matrix can be written as ρt The state |φ t 〉 can be represented as a matrix product state (MPS) [95], and the single-site operator e iθ ŝx j can be expressed as a two-by-two matrix.Utilizing this representation, the moment generating function can be computed by sandwiching the operators between the matrix product states, as depicted in Figure 12.By obtaining G l (θ , t), the complete probability distribution P l (m, t) is computed numerically by discretizing the Fourier integral in equation 8.
Figure 12: Computing the generating function G l (θ , t) in matrix product state representation.The site i is chosen such that the subsystem of size l is in the center of the full system.We conducted two types of error analysis to assess the accuracy of the numerical results.In figure 13, we assess the absolute error of the TDVP algorithm in comparison with the numerically exact full diagonalization results for a system with size N = 14 and various post-quench parameters.Figure 13, panel (a), shows the absolute error in the energy density of the thermal states, defined as

B.3 Errors and data convergence
The absolute error remains of the order O(10 −5 ) or smaller across the entire temperature range under consideration.Figures 13, (b), (c), and (d) show the absolute error in domain wall kinks, defined as |〈 k〉 ED − 〈 k〉 TDVP |, following a quantum quench to various post-quench parameters.〈 k〉 is defined in equation 13 and the computational basis {|n〉} is its simultaneous eigenbasis.Notably, the error rapidly converge and is of order O(10 −6 ) or smaller for all the cases studied.
In Figure 14, we investigate the convergence of the TDVP data for DT(t) by computing the relative errors |DT χ 1 (t) − DT χ 2 (t)| for three increasing bond dimensions.Our observations reveal that the relative error eventually converges and consistently remains in the order O(10 −3 ) or smaller for all cases.It is noteworthy that the error for α = 0.0 is several orders of magnitude smaller than that for the other values of α.This is attributed to α = 0.0 being an integrable point with an extensive number of conserved quantities, and therefore has a smaller Hilbert space to be explored compared to non-integrable points.Furthermore, for α = {1.5, 1.9}, the error for h f = 0.3 is approximately two orders of magnitude smaller than that for h f = 0.6.This discrepancy arises because the former case exhibits dynamical confinement, which effectively suppresses the spread of correlations and constrains the total Hilbert space that can be explored during time evolution.For values of α > 2, the long-range Ising model falls within the regime of short-range interactions and does not exhibit any finite-temperature phase transitions [89].Extensive investigations into the critical properties of the thermal phase transition in the quantum longrange Ising model have been conducted using numerically exact path integral Monte Carlo methods [90].The thermal phase transition is qualitatively depicted in Figures 15 for specific parameter values: α = 1.5, 1.9 and h = 0.3, 0.6.As described in Section B.1, the simulation begins with a maximally mixed state at β = 0.This initial state is characterized by a sharply peaked Gaussian distribution of P(m) centered around m = 0, which signifies a strongly para-magnetic phase.As the system is gradually cooled by increasing β, the distribution gradually widens, eventually becoming nearly flat around the critical temperature.A further reduction in temperature leads to the emergence of a bimodal distribution of P(m), which is indicative of the ferromagnetic phase.Notably, this transition from a unimodal Gaussian distribution to a bimodal distribution highlights the 2 symmetry that is inherent in the long-range Ising Hamiltonian.Confinement phenomena in the long-range Ising model result from ferromagnetic interactions extending over long distances between the interacting spins.However, the strength of confinement varies within different regions of the phase space [44,45].In this section, we present comprehensive numerical results pertaining to the temporal spreading of correlations in the long-range Ising chain following a sudden quench to various post-quench Hamiltonians starting from a fully polarized initial state denoted as |ψ i 〉 = |←, ←, . . ., ←, . . ., ←, ←〉 x .In panels (a), (b), and (c), we examine a fixed value of α = 1.9 while varying the transverse field h = 0.3, 0.6, 0.8.Notably, panel (a) shows a pronounced signature of confinement, which gradually diminishes as the value of h increases, as shown in panels (b) and (c).This behavior is expected because the transverse field competes with long-range interactions and weakens the confinement effect.In panel (d), we observe a linear light cone spreading of the correlation with the maximum possible velocity, v max = 2h [42].

Figure 1 :
Figure 1: (a) Global Quench protocol: System is initialized as the ground state of a trivial hamiltonian Ĥi (in our case the initial state is a Greenberger-Horne-Zeilinger (GHZ) state), at time t = 0 the system is suddenly quenched to a final Hamiltonian Ĥ f and the initial state is unitarily evolved with the final Hamiltonian.(b) The nonequilibrium state following a global quantum quench can exhibit different relaxation behavior; Path 1: a direct relaxation to thermal equilibrium with a single time scale, Path 2: a quick relaxation to a long lived prethermal state eventually followed a relaxation to thermal equilibrium, Path 3: a strong retention of initial memory and suppression of relaxation to thermal equilibrium.

Figure 2 :
Figure 2: Numerical extraction of β eff corresponding to a global quantum quench.The horizontal black dashed line represent the energy density attributed to the quench.The colored lines represents the energy density as the function of inverse temperature β for the corresponding post-quench parameter (in legend).The point at which the colored lines intersects the black dashed lines represents β eff for the corresponding post-quench parameters (represented by vertical colored lines).

Figure 5 :
Figure 5: Time evolution of domain wall kinks, 〈 k〉, following a global quantum quench to three interaction strength values α ∈ {0.0, 1.5, 1.9} (Panels (a), (b), and (c) respectively) and h f = 0.3 and h f = 0.6.The dashed horizontal lines represent the expected thermal value of domain wall kinks, 〈 k〉 TH corresponding to the quenches.

Figure 6 :
Figure 6: Maximally mixed density operator at β = 0 as the tensor product of local identities.

Figure 8 :
Figure 8: One half of the ρβ=0 in the locally purified form.

Figure 11 :
Figure 11: Expectation of the local operator Ôi in thermal density operator ρβ .

Figure 13 :
Figure 13: Absolute error in the energy density,|ε ED β − ε TDVP β |, of thermal states -(a).Absolute errors in the evolution domain wall kinks |〈 k〉 ED − 〈 k〉 TDVP | following a quantum quench -(b),(c),(d).The numerically exact results are calculated using equations 19 and 20 as detailed in A. TDVP results are obtained with bond dimension χ = 128.The system size considered is N=14.

Figure 14 :
Figure 14: Convergence of the TDVP data for DT(t) with increasing bond dimensions, χ = 60, 90, 128, for six different post-quench parameters considered in the main text.The black dashed line is for visual guidance.

Figure 15 :
Figure 15: Thermal phase transition of long range Ising model at four different points in parameter space.The initial state in all cases is the maximally mixed state at infinite temperature represented by ρβ=0 , refer to 6.The color coding from red to blue signifies decreasing temperature.

Figure 16 :
Figure 16: Real time dynamics of half chain connected correlation function ŝx k ŝx k+∆ c after a global quantum quench of the transverse field starting from a fully polarized initial state.The dashed black lines is v max = 2h line for nearest neighbor transverse field Ising model [42].

Figure 16
Figure 16 illustrates the time evolution of the half chain connected correlation function ŝx k ŝx k+∆ c = ŝx k ŝx k+∆ − ŝx k ŝx k+∆ in a chain of 200 spins, where k is kept fixed at the center of the chain.In panels (a), (b), and (c), we examine a fixed value of α = 1.9 while varying the transverse field h = 0.3, 0.6, 0.8.Notably, panel (a) shows a pronounced signature of confinement, which gradually diminishes as the value of h increases, as shown in panels (b) and (c).This behavior is expected because the transverse field competes with long-range interactions and weakens the confinement effect.In panel (d), we observe a linear light cone spreading of the correlation with the maximum possible velocity, v max = 2h[42].