Decoherence and recurrences in a large ferromagnetic F = 1 spinor Bose-Einstein condensate

Decoherence with recurrences appear in the dynamics of the one-body density matrix of an F = 1 spinor Bose-Einstein condensate, initially prepared in coherent states, in the presence of an external uniform magnetic field and within the single mode approximation. The phenomenon emerges as a many-body effect of the interplay of the quadratic Zeeman effect, that breaks the rotational symmetry, and the spin-spin interactions. By performing full quantum diagonalizations very accurate time evolution of large condensates are analyzed, leading to heuristic analytic expressions for the time dependence of the density matrix, in the weak and strong interacting regimes. We are able to find accurate analytical expressions for both the decoherence and the recurrence times, in terms of the number of atoms and strength parameters, that show remarkable differences depending on the strength of the spin-spin interactions. We discuss the nature of these limits in the light of the thermodynamic limit.


Introduction
The observation of stationarity in quantum systems relies on the existence of pair wise collisions in large conglomerates of atoms, either in an isolated environment or in contact with a larger environment [1][2][3][4][5][6][7][8][9][10]. While in the former case the process through which the system reaches the stationary state may be termed intrinsic decoherence [11,12], in the later case it is simply known as decoherence. Although measurements of purity have been used to search for quantum coherence loss [13], the most common scheme to study this phenomenon is by tracking the time evolution of the density matrix describing the state of the body. In the case of isolated systems composed of a macroscopic number of constituents N, the reduced density matrices ρ (n) with n = 1, 2, ... < N, are the proper ones to record for signatures of stationarity and in contrast, in systems with an equal number of constituents N but embedded in larger environments of size N N, the corresponding quantity is the N−body density matrix ρ (N) . In both cases, the decoherence time is the interval elapsed until the system reaches the stationary state, that is, the time in which the average value of the physical quantities do not vary any more. Decoherence, intrinsic or not, is a many-body phenomenon appearing in few-body or small parts properties of a system, and should depend on N. Recurrences typically appear both because of the finiteness of the system and of the associated Hilbert space. This type of decoherence and eventual relaxation to thermal equilibrium is what occurs in a closed macroscopic body [14].
Decoherence and recurrence phenomena have been widely studied both, experimentally and theoretically. Within the experimental side some examples are the dephasing in interference fringes measured in condensates, the decay of laser induced polarization in spectroscopic experiments [15,16], amplitude damping in qudit photonic states [17], and decoherence induced in single molecule junctions [18] among others. Regarding the theoretical description, dissipation or relaxation processes are usually included by means of a master equation that take into account the coupling of the atoms with the environment. The literature in this respect has a very long history and it is extremely abundant, see Refs. [1][2][3][4][5][6][7][8][9][10] as notorious examples.
In this article, within a full quantum scheme, we explore the phenomenon of intrinsic decoherence as well as the appearance of recurrences, in the time dynamics of an spinor F = 1 Bose-Einstein condensate (SBEC) composed of a mixture of three different hyperfine spin components, in the presence of a uniform magnetic field. Depending on the sign of the spin-mixing interaction strength, the atomic cloud can be polar if positive, such as a gas of 23 Na atoms, or ferromagnetic if negative, as in a gas of 87 Rb atoms [19]. We shall study here only the ferromagnetic case since the condensate acquires a macroscopic spin texture that makes it behave as a behaves as a "giant" spin. It is worth mentioning that the system we address is very similar to the recent experimental investigation on the spin dynamics of an F = 1 87 Rb spinor macroscopic condensate [20]. The present model has also been used to study quantum phase transitions in space [21]. The dynamics of the spin mixture is followed in the presence of the external magnetic field that gives rise to linear and quadratic Zeeman contributions, with strengths p and q, which together with the interaction term, of strength η, bring the system towards stationarity for any arbitrary initial condition.
Since we are able to very accurately diagonalize the Hamiltonian of the system up to N ∼ 10 4 , we can study the dynamics of any initial state and then calculate the reduced density matrix. Here, we analyze the one-body density matrix for the full family of coherent states [22], as one expect them to be the most resilient to the spin interaction. And indeed, if there is no magnetic field, the coherent states show no decoherence, but as soon as the full rotational symmetry is broken by the presence of the quadratic Zeeman interaction, decoherence and its associated recurrences ensue. After a long set of calculations we were able to find out what amounts to be the main contribution of this article. On the one hand, we found out that there is a dimensionless parameter N|η|/q that clearly separates three regimes, a weak interacting one N|η|/q 1, a strong one N|η|/q 1 and a crossover N|η|/q ∼ 1, in which, although all show decoherences, stationary states and recurrences, their dependence on the parameters q and η and specially on N, are very different. But more importantly, after the thorough and lengthy numerical study, and with the insight of the analytic solution to the non-linear single mode model [23], we were able to heuristically deduce analytic expressions for the full one-body density matrix in the weak and strong interacting regimes. These expressions are certainly the leading order contributions of the full unknown analytical solution and show a remarkable agreement with the predictions of the full quantum numerical solution. We are thus able to provide analytical expressions for the recurrence and decoherence times in terms of the parameters q, η, N and their dependence on the particular initial states.
This manuscript is organized in 5 sections. In section 2 we introduce the model Hamiltonian that describes the spinor BEC within the single mode approximation (SMA) and discuss how we are able to accurately obtain its full quantum diagonalization. In section 3 we introduce the one-body density matrix its time evolution for the full family of coherent states in the ferromagnetic case, preparing the stage for section 4 that shows our main contributions regarding the discussion of decoherence recurrences in the weak and strong interacting regimes. Finally a discussion and a summary of this work is presented in section 5.
As mentioned in the Introduction, we consider an F = 1 spinor Bose-Einstein condensate (SBEC) in the presence of an external uniform magnetic field in theẑ direction, breaking the full rotational symmetry of the atomic cloud. Under these conditions it is valid to assume that all spin states have the same spatial wave function yielding a decoupling with the spin degrees of freedom, allowing for a study of the dynamics of the spin modes only. This is the single-mode approximation (SMA) [19,[24][25][26][27][28][29]. Hence, while one looses the spatial dependence of the condensate, one can perform a many-body full quantum analysis. Although not shown here for brevity, we have performed full evolution of an F = 1 Gross-Pitaevskii equation, in the ferromagnetic regime, and found that SMA holds perfectly well for uniform magnetic fields, see Ref. [30] for details of our method.
Implementing the above approximation yields the SMA Hamiltonian of an F = 1 SBEC, [24,25,31] The first and second terms represent the linear and quadratic Zeeman contributions, while the third one is the spin-mixing interaction. The above operators are, in turn, given in terms of the creation and annihilation operators of bosonic atoms in the spin states +1, 0 and −1, in obvious notation, Certainlyˆ f = (f x ,f y ,f z ) is the many-body angular momentum vector operator, witĥ The coupling parameters p, q and η denote the strengths of the respective contributions which, in principle, can be experimentally tuned. As already mentioned, η can be either positive or negative, indicating the polar and ferromagnetic nature of the interactions [19]. The dynamics that we study below will be limited to the ferromagnetic case, η < 0. Introducing the spin state number operatorsn σ =â † σâσ , one can also writef z =n +1 −n −1 andQ =n +1 +n −1 , forms that can be useful in interpreting our results below.
For a given number of atoms N the size of the Hilbert space is Ω = (N + 1)(N + 2)/2 and, therefore, the size of Hamiltonian given by equation (1) scales as ∼ N 2 × N 2 . In order to find the time evolution of the system we have to diagonalize the Hamiltonian, a difficult task that can be eased by exploiting its symmetries. Obviously, the total number operatorN =n +1 +n 0 +n −1 commutes with the HamiltonianĤ. In addition, due to its Lie structure, it is easy to show that f z ,Ĥ = 0. Hence, instead of using the "natural" basis of number occupation |n +1 , n 0 , n −1 , in an obvious notation, one finds a better alternative to use |M, n 0 , with M = n +1 − n −1 the eigenvalues off z , with values −N, −N + 1, . . . , N − 1, N. Atom number conservation N = n +1 + n 0 + n −1 yields the third quantum number, obviated in the state labels. In this basis, the Hamiltonian is block diagonal in M and the matrix elements show simple expressions, where n +1 = (N + M − n 0 )/2 and n −1 = (N − M − n 0 )/2. We note that the matrix elements in Eqs. (5)-(7) do not link blocks with different values of M, however, they connect different matrix elements with jumps of 2 for n 0 . Then, the Hamiltonian H can be diagonalized and written as, where, for each block M, m M = 1, 2, . . . , m max M , and m max M = (N − |M| + 1)/2 or m max M = (N − |M| + 2)/2 if M is odd or even, respectively. This analysis can also be applied to SBEC with F > 1 [32,33]. The block-diagonal structure of the Hamiltonian given by Eq. (1) allows not only to very accurately calculate the energy eigenvectors and eigenvalues for values up to N ∼ 10 4 , as given by Eq. (8), but also to implement full quantum evolution of any initial quantum state for arbitrary values of time. For this, we perform numerically exact diagonalization block by block using pyCUDA linear algebra routines, in a server with 8 CPU Core, 128 GB RAM and with a graphic card Nvidia Tesla C2075. Since this diagonalization can be obtained for a wide variety of parameters, the ensuing time evolutions we obtain do not suffer from accumulation of errors, thus maintaining the same numerical precision at any time step. Depending on the value of N and the needed time steps, the full evolution of a given initial state elapsed between few GPU seconds up to several hours. As an example easy to visualize, in Fig. 1 we display the Hamiltonian structure for N = 6 particles with a Hilbert space size of 28, each blue square representing a block of magnetization M, whose size is given below. The intensity of color blue and the size of the blocks depends on the value of M. Although we do not address it here, we show in Fig. 2 the spectra and its degeneracy for different values of p, q and η for N = 10 3 , to illustrate its richness.  Despite the fact thatf z andQ commute, the interesting and rich behavior of this model arises because of the non-conmutativity of the quadratic Zeeman term ∼Q and the spin interaction ∼f 2 . Further, in the absence of the quadratic Zeeman interaction, q = 0, the Hamiltonian can be analytically diagonalized both for F = 1 [34] and F = 2 [32]. The presence of the quadratic term, q = 0, breaks the axisymmetry [31] and this system becomes an excellent one to study Quantum Phase Transition [35], quench-dynamical behaviors [36], Quantum Kibble-Zurek Mechanism [37], spin fragmentation [38] among other mechanisms. On the contrary, if there is no atomic interactions, η = 0, the problem is also trivial and right away diagonal as given by Eqs. (5)-(7). Hence, reiterating, the presence of the atomic interactions is mediated by the presence of the Zeeman quadratic term.

Time evolution of the one-body density matrix in coherent states
With the full quantum diagonalization described above we can calculate the time evolution of arbitrary initial states |Ψ 0 . Here, we will concentrate in the full family of initial coherent states. For this purpose we briefly mention that, as expected, the time evolution can also be performed by blocks in the following manner. First, since the unitary propagator operator can be written as, one just needs the overlaps M, m M |Ψ 0 to find the time evolution of any |Ψ 0 . However, since many of the usual physical quantities are typically one-or two-body operators, we shall devote our attention to the one-body density matrix, which allows for calculating all the statistical properties of all one-body operators, such as the matricesf x ,f y , andf z . For this we note that since any one-body operator can be written as,Ô with k and j taking the values +1, 0 and −1 and O ij complex numbers, the expectation value of any operatorÔ (1) requires the knowledge of the expectation values of the operators a † k a j for all values of kj. These expectation values are those of the one-body density matrix. Explicitly, for a given initial state |Ψ 0 , the one-body density matrix for all times is given by, Although we do not exploit here, it is worth mentioning that if we limit ourselves to one-body properties, one could also use the properties of the Gell-Mann spin 1 matrices [39].
As we can accurately calculate the matrix ρ kj (t) for any time using Eq.(11), we now turn our attention to the initial set of coherent states, [22] where ζ 1 = e −iϕ cos θ , with θ and ϕ the usual angles of the unit sphere, and |0 denoting the vacuum state with no particles. Alternatively, a coherent state can also be written as, A very important property to take into account is that these coherent states are eigenstates off 2 , that is, f 2 |θ, ϕ = N(N − 1)|θ, ϕ . The main features of the time evolution of these states, as we amply discuss and show below, is that the elements of the one-body density matrix in these states show an initial Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 12 November 2020 doi:10.20944/preprints202011.0353.v1 oscillation that suffers intrinsic decoherence followed by a stationary state and revivals at later times, with this behavior being repeated ad infinitum. It is evident that both the decoherence and recurrence times depend on the p, q and η parameters as well as on the initial state (θ, φ), but an important issue is its dependence on N. As mentioned in the Introduction we have indentified that the dependence on N is very different in two opposite limits N|η|/q 1 and q/N|η| 1, evidently called weak and strong interacting regimes.
Since one can show that the set of operatorf z ,Q andf 2 are not part of a Lie algebra, the finding of an analytic expression for the reduced density matrix appears as a very difficult task. Nonetheless, the main contribution of this article, is to show that the time evolution of the SBEC one-body properties can be summarized quite precisely with explicit (semi) analytic expressions for the time evolution of the density matrix elements ρ kj (t), given in Eq. (11), for arbitrary initial coherent state such as Eq. (12), in the weak and strong limits. These expressions, as well as their validity limits, are found in a heuristic manner based on a very large number of precise numerical evaluations of time evolutions for a wide variety of values of the Hamiltonian parameters and for a collection of different initial coherent states.
In order to introduce our expressions for the one-body density matrix in the following section, we first discuss preliminary exact results. Note that the one-body density matrix ρ kj (t), given by (11), can be first expressed as, where and |θ ≡ |θ, 0 is the coherent state for any value of θ and ϕ = 0. We see that for any initial state the role of the pf z term and the angle ϕ simply amount to a phase factor, while the whole dynamics is ruled by the quadratic Zeeman qQ and the interaction ηf 2 terms. Whether one term or the other dominates depends on the relative values of the strength parameters q and η. We further note that if q = 0, the interactions play no role since the coherent states are eigenstates off 2 . Therefore, in order to observe the effect of both terms, both η and q must be nonzero.
As an illustration of typical time evolution of coherent states, in Fig. 3 we show sequences of decoherences and recurrences in the weak N|η|/q 1, strong q/N|η| 1 and crossover N|η|/q ∼ 1regimes. We have found that in the extreme cases the recurrences appear at periodic intervals, thus making feasible their prediction, while in the crossover, as the two effects contribute similarly, the sequences can be very irregular and we make no attempt to analyze them here.

Decoherence and recurrences in the weak and strong interacting regimes
To proceed with the finding of the weak and strong interacting extremes, we note that the density matrix, Eq. (14), can be written in two alternative forms, ρ kj (t) = θ|e − ī h η τ 0V I (τ)dτ a † k a j e ī h η τ 0V I (τ)dτ |θ e i(k−j)pt/h e iδ k,j±1 (δ j,2 −δ k,2 )qt/h (16) or whereV and e T stands for the time-ordered exponential. The first form, Eq. (16), is appropriate for the study of the weak interaction N|η|/q 1, while the second one, Eq. (17), for the strong one q/N|η| 1. Note that the exact limits η = 0 and q = 0 are explicitly recovered in Eqs. (16) and (17).
In addition to the analysis of a large number of numerical calculations of the time evolution of coherent states, varying all the relevant parameters, we also used as insight the analytical solution of the non-linear single-mode Hamiltonianĥ = −µâ †â + uâ †â †ââ , see Refs. [40] and [23]. From these is not difficult to conclude that the envelope of the decoherences, at all recurrence times, is of a gaussian form. Furthermore, the recurrences are observed to appear at periodic times. We proceed now to present the heuristic (semi) analytic forms of the density matrix in the two extremes. Figure 4 shows typical behavior of elements of the one-body density matrix in the weak interaction regime, for a particular case, see figure caption. We find that all essential features of the time evolution of the one-body reduced density matrix can be fitted quite well by the following explicit expressions.  ρ jk (t) ≈ ρ jk (0)F jk (N, η, θ, t)e −2i(j−k)ηN cos θt e −i(j−k)pt e iδ j,k±1 (δ k,0 −δ j,0 )qt j = k where F jk (N, η, θ, t) = e g jk (θ)N(cos(2 f jk ηt)−1) for all j, k .

Weak interaction regime N|η|/q 1
Above, j, k = −1, 0, +1; f ij are symmetrical with f 1,0 = 1, f −1,1 = 2, f −1,0 = 3 and all diagonal equal, f jj = 4. We were not able to find analytic expressions for the coefficients g jk (θ) and A jj (θ) but they can be numerically fitted, as we show them in Figs. 5 and 6. The diagonal terms g jj (θ) are all equal, shown in Fig. 6. The initial condition is given by the matrix: which was obtained using the lie algebra method [41] •  Figs. 4 and 7 show the density matrix ρ jk (t) and the expectation value of the operatorf x for a single example in the limit N|η|/q 1 (η = −10 −5 and q = 7 in dimensionless units, for N = 300). These figures illustrate the excellent agreement of the heuristic expression given by Eq. (19) with the accurate numerical calculations. Fig. 4(a) shows both cases for the initial decoherence, while Figs. 4(b), (c) and (d) the agreement in predicting the occurrence of the recurrences, for ρ −1,0 , ρ −1,1 and ρ 0,0 . Fig. 7 shows the behavior of the elements (a) ρ 1,0 and (b) ρ −1,0 . As can be seen in Eq. (19), these two elements oscillate with a very large frequency q = 7 and with a very small correction ±2Nη cos θ, each taking one of the signs. This tiny difference between the oscillations of ρ 1,0 and (b) ρ −1,0 shows itself in the beating pattern of the expectation value off x (see Eq.(4)) shown in Fig. 7(c) and (d): the former is the exact numerical calculation and the latter the heuristic expression given by Eq. (19). We found a quite remarkable agreement for any value of θ.
The behavior of intrinsic decoherences, followed by quasi-stationary states to further revivals or recurrences in a periodic fashion, is essentially contained in the function F jk (N, η, θ, t) ≈ Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 12 November 2020 doi:10.20944/preprints202011.0353.v1 . The beating pattern in f x (t) is due to a very small difference of the high frequency oscillations of oscillations between ρ +1,0 and ρ −1,0 , see Eq. (19). The (dimensionless) parameters are N = 300, η = −10 −5 ; q = 7, θ = 7π/30, p = 0 and ϕ = 0. e g jk (θ)N(cos(2 f jk ηt)−1) within the density matrix, see Eq. (19). The function F jk indicates right aways that there are recurrences at periodic intervals, τ n = nτ rec with n = 1, 2, 3, . . . , and that we label as the weak-interacting recurrence time; it does depend on the particular matrix element but the main point is its inverse proportionality to η and its independence on N. Then, at each recurrence, starting with the initial time t = 0, the evolution appears to decohere in a shorter time scale. This can be found by approximating the exponential for short times, thus identifying the weak-interacting decoherence time, e −t 2 /2τ 2 dec , In this case, in addition to the dependence on η, there appears a dependence on the number N of atoms. We defer further discussion of these time scales to the last section.

Strong interaction regime q/N|η| 1
In this regime, contrasting with the weak one, the leading terms depend mostly on the quadratic Zeeman coupling q and the one-body reduced density matrix can now be very precisely fitted by a much simpler form, with ρ kj (0) given by Eq. (22). Once more, the structure of periodic recurrences with decoherences in relatively shorter times is given by the real exponential term exp( N 2 sin 2 θ cos q(j−k) hN t − 1 ), indicating that at the times τ n = nτ (q) rec , n = 0, 1, 2, . . . , periodic recurrences occur, We label this as the strong-interaction recurrence time. The strong-interaction decoherence time can be readily find by expanding the previous exponential at short times, yielding again a gaussian function Note that both times τ (q) rec and τ (q) dec grow as does the number N of atoms. We also defer further discussion of this behavior to the next section.
The diagonal terms ρ jj (t), for j = −1, 0, +1, are essentially constant equal to their initial values at t = 0. The precision of our calculations allows to set, ρ jj (t) ≈ ρ jj (0) + a jj (θ)e with a jj (θ) a term that it is at most of the order of 10 −9 , compared with the initial value ρ kj (0), for In Fig. 8 we show an example of the evolution of the density matrix in the strong interacting regime q/N|η| 1 ( q = 1, p = 1 and η = −30000 in dimensionless units for N = 700). This illustrates both the typical behavior of decoherences and recurrences and the agreement with the heuristic fitting given by Eq. (26). In panels (a), (b) and (c) we show the evolution of the real part of the off-diagonal elements of the density matrix that show the sequence of decoherences followed by stationary states and recurrences. In panel (d) we show the agreement with the predicted oscillations in terms of the q and p parameters.
As shown in Fig. 9, and confirmed by Eq. (26), the behavior of one-body observables, such as f x , show a very detailed and predicted structure of decoherences and recurrences, with the same stationary state for a very large parameter space. As seen in this and in all the above figures and expressions, the one-body density matrix of the stationary state is diagonal in the spin basis (−1, 0, +1) along the z−direction. Note from the figure that if q 1, but η finite, the recurrences can be separated by very long times, that increase as N grows, hence, yielding a true stationary state in the thermodynamic. This is further discussed below.

Discussion and Final Remarks
The agreement between the proposed heuristic expressions for the time evolution of the one-body density matrix compared with full (numerically) exact full quantum calculation is certainly quite good, allowing us to validate our conclusion that those expressions are the leading non-trivial order, becoming better as N grows. Once more, in both cases the general behavior shows periodic (intrinsic) decoherences of an oscillatory behavior, into stationary states, with periodic recurrences. However, the precise nature of this phenomenon is very different in the two limits considered. The most notorious difference between the extreme regimes resides on the fact that the role of the number of particles N is very different in each case. Note that in the weak side τ rec ∼ N. It is of interest to mention that the non-linear single-mode approximation, as described in Refs. [40] and [23], corresponds to the weak-interaction regime. From a priori perspective of the behavior of a macroscopic body, N 1, one can state that the strong-interaction limit indeed fulfills such an expectation. That is, from a thermodynamic point of view, a body should take longer to reach a stationary, or equilibrium state, the larger it is. Hence, τ (q) dec ∼ √ N appears as a reasonable result. At the same time, one expects that the recurrences should take longer to appear as the body becomes larger and, again, τ (q) rec ∼ N appears as a good conclusion. It is the more appealing the fact that the ratio τ dec /τ rec ∼ 1/ √ N, a typical condition for thermodynamic stability of a macroscopic system [14]. That is, although both times diverge in the thermodynamic limit, in the appropriate scale, the decoherence time tends to zero as ∼ 1/ √ N compared to the recurrence one. Therefore, we find very interesting to observe that while the decoherence and recurrence times in the weak-interaction case, N|η|/q 1, do not follow the "expected" thermodynamic trend, still we observe that the ratio τ (η) dec /τ (η) rec ∼ 1/ √ N scales appropriately. In any case however, although the role of the number of atoms N is very different in each limit, still the ultimate responsible for the decoherence and recurrence phenomena is the interaction among the constituents of the body. We believe this strong difference should also be present in the structure of the stationary states attained between consecutive recurrences; that is, we expect states closer to typical thermal equilibrium in the strong rather than in the weak interaction regime. This deserves a separate and detailed study. In the crossover regime N|η| ∼ q both terms compete, their effect is intertwined and while the system indeed shows decoherences and recurrences, the latter no longer occur at prescribed times depending on either √ N or N, see Fig. 3 (c). To summarize, we first highlight that one can exactly (numerically) diagonalize the Hamiltonian of an interacting F = 1 SBEC in the SMA approximation for large number of atoms and that the method can be extended to larger spins F > 1. This allows us to probe the full quantum dynamics of any initial state. For the purposes of the present study, we have chosen here the relevant family of the coherent states. We have studied the corresponding reduced one-body density matrix and have found heuristic analytical expressions that fit remarkably well its dynamics in the weak N|η|/q 1 and strong N|η|/q 1 regimes, indicating the interplay among the non-linear Zeeman effect ∼ q and the pairwise spin interaction ∼ η, mediated by many-body effects ∼ N. Our expressions predict the decoherence time and the recurrence periods in terms of these quantities and become more precise as N grows. Since the corresponding unitary propagator cannot be analytically found, due to its lack of closed Lie algebra, our results may indicate a path to find it in a series whose leading order term agrees with the heuristic expressions here found.
There are many additional venues to further explore. We mention two here. One is the study of two-body correlations. For this we need to calculate the reduced two-body density matrix, which with our method does not appear as a very difficult task. It would be very interesting to find how these correlations behave along the different regimes. In addition, a second and complementary study is the elucidation of the nature of the stationary state between recurrences. In particular, one could enquire if there is an "equilibration" similar to the Eigenstate Thermalization Hypothesis (ETH) [42,43]. That is, since we can also consider thermal states with a Gibbs distribution [14], due to our exact diagonalization, we could compare expectation values of physical quantities, such as angular momentum and their correlations, obtained from thermal states with those obtained in single states in the stationary regions of their evolution.