Intrinsic 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 one-body density matrix, in the weak and strong interacting regimes, for initial coherent states. 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. The features of the stationary states in both regimes is also investigated. We discuss the nature of these limits in the light of the thermodynamic limit.


I. 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. In reality, since one can consider the system (S) under study and the environment (R) as a composite closed system (S + R), their time evolution is always unitary and the observation of decoherence always refers to reduced quantities, namely, observables of much fewer degrees of freedom of the total ones of the isolated system, say N S N R , with N the number of degrees of freedom. Therefore, when studying a system in contact with a reservoir, one deals with the reduced density matrix of the former, ρ S , integrating out the environment. This leads typically to master equations for the reduced density matrix of the system that has already ingrained the decoherence effects induced by the interaction with the environment [2,5,6,[8][9][10]. In general, 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 [13,14], the amplitude damping in qudit photonic states [15], and the decoherence induced in single molecule junctions [16] among others. Measurements of purity have also been used to search for quantum coherence loss [17].
On the other hand, if the system under study is a large one, N 1, but isolated from any environment, even though the evolution is always unitary, the system does relax or decohere to a stationary state [18] within a quantifiable decoherence time, in the sense that expectation values of one-body observables, such as temperature, magnetization and, say, few-body correlations, behave as if the full system had relaxed to a stationary state for times longer than the decoherence time. If the system is finite, there appear recurrences or revivals of states near the initial one, but in typical systems the recurrence time may grow with no bound in the thermodynamic limit N → ∞. The behavior of those few body properties can be directly studied solely with their corresponding reduced few-body density matrices, whose time evolution is no longer unitary within their reduced Hilbert space. It is of interest for our purposes that the experiments in ultracold gases showing Bose-Einstein condensation [19][20][21][22] are the closest in dealing with true isolated systems. Yet, these gases, due to atomic collisions, do relax to equilibrium and, when in the presence of external magnetic fields do show decoherence phenomena [23,24]. The study of the magnetization of the latter case is the subject of the present article. Hence, to be specific, we call intrinsic decoherence to that of few-body properties, in an otherwise very large isolated system. Recurrences are not observed in the ultracold gases since unwanted processes, such as three-body collisions, make the thermal state unstable in a relatively short time [25]. An interesting question concerns the nature of the stationary state and the enquiry if it is a thermal state or not. In this regard the essence of the Eigenstate Thermal Hypothesis (ETH) [26][27][28][29][30] is to establish that the thermodynamic properties of few body observables are contained in the eigenstate closest to the equilibrium average energy. Hence, a simple test is to compare the few-body density matrices of the stationary state of the actual unitarily evolving state with those of the energy eigenstate with an energy similar to the mean of the state evolving in time.
Deviations from this typical behavior, such as the many-body localization phenomenon [31], are also of current interest. In the light of these observations, we advanced here that depending on the two-body interaction strength in an F = 1 spinor Bose-Einstein (SBEC) condensate in the presence of external homogeneous magnetic fields, we do observe intrinsic decoherence of the magnetization but with different features leading to both, typical and non-typical stationary states.
To be precise, within a full quantum scheme, we study here 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, starting in a well defined coherent state. 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 [32,33]. We shall study here the ferromagnetic case only since the condensate acquires a macroscopic spin texture that makes it behave 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, where use of the SBEC as a sensible magnetometer is explored [24]. The present model has also been used to study quantum phase transitions in space [34]. 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 η, give rise to the phenomena here analyzed.
Since we are able to very accurately diagonalize the Hamiltonian of the system up to N ∼ 10 4 atoms, we can study the dynamics of any initial state and then calculate reduced density matrices of few bodies. Here, we analyze the one-body density matrix for a full family of coherent states, as one expects them to be the most resilient to the spin interaction, clearly showing Larmor-like oscillations in the expectation value of the measurable spin (or magnetization) one-body observable. We point out that coherent states are usually expected to yield the closest to quasi-classical dynamics of the magnetization, in a mean-field fashion, with little or no decoherence [35]. And indeed, if there is no quadratic Zeeman coupling, the coherent states show no signs of decoherence, but as soon as the full rotational symmetry is broken by the presence of such a coupling, this, in conjunction with the natural two-body collisions yield decoherence or collapse of the very definite Larmor oscillations generated by the presence of the external field. Due to the finite number of atoms and the relatively small Hilbert space size, the oscillations recur or revive after longer periods of time, to decohere again. After a long set of calculations we were able to find out three clear regimes depending on the value of the dimensionless parameter N |η|/q, 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 and recurrences, their dependence on the parameters q and η and specially on N , as well as the nature of their corresponding stationary states, are very different. Also very importantly, after the mentioned numerical study, and with the insight of the analytic solution to the non-linear single mode model [36,37], 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. We also highlight the fact that the strong interaction regime shows a "typical" behavior of a macroscopic system since the decoherence and recurrence times show a expected dependence on N , in the sense that as N → ∞, the recurrences time grows with no bound, thus making the quasistationary state closer to a true one. In addition, the nature of the stationary state behaves similarly to ETH as having the same reduced one-body density matrix as the eigenstate whose energy equals the gas average energy.
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 and its time evolution for a 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.
The many-body F = 1 SBEC Hamiltonian with linear and quadratic Zeeman couplings to an external homogeneous magnetic field B, within the contact approximation, is whereψ α (r) are annihilation operators of particles at r with spin α = −1, 0, +1 and F αβ are the F = 1 angular momentum matrices. c 0 and c 2 are interaction coefficients proportional to the corresponding s-wave scattering lengths. If c 0 > c 2 the system is polar and for c 0 > c 2 ferromagnetic [32]. U (r) is an external confining potential, typically harmonic. In general, the field operator is given byψ where φ mα ( r) are elements of a basis of the one-particle Hilbert space andb mα the corresponding creation operators. For ultracold gases, a usual approximation is to consider a self-consistent to be determined ground state wavefunction only Ψ 0α (r, t), such that ψ α (r, t) ≈ Ψ 0α (r, t)b 0 , that is all spin states are in the same ground state. This leads to a set of three coupled Gross-Pitaevskii (GP) equations [32,33]. These can be solved numerically and, as shown in Fig. 1, the solution for a homogenous magnetic field shows that the spatial part is unaffected by the presence of such a field, showing Larmor-like oscillations of the magnetization, see below. That is, the dynamics occurs only in the spin degrees of freedom and it is not transferred to spatial excitations such as phonons and vortices. This would not be so if the external field is inhomogeneous [38]. The previous observations indicates that a single-mode approximation can be made at the level of GP description (SMA-GP), that is, an assumption that the ground state wavefunction is the same for all spin components α .
So, if we set Ψ 0α (r, t) = ψ 0 ( r)ζ α (t) with ζ α (t) a simple 3-component spinor, we can integrate out the spatial part and find very simple dynamical equations, called SMA-GP (not shown here), for the time evolution of ζ(t). The result of those is shown also in Fig. 1 where we plot the time evolution of the magnetization f (t) = drΨ * 0α (r, t)F αβ Ψ * 0α (r, t) with full 3D GP, see Ref. [39] for details of our methods, and f (t) = ζ * α (t)F αβ ζ * α (t) with the SMA-GP equations. The agreement is essentially perfect. The unphysical feature of the SMA-GP equations, as well as of the full 3D GP, is that they are incapable of showing decoherence effects. This is not surprising since GP equations assume that there is a single macroscopic wavefunction for the ground state even in the presence of dynamical effects. with a field B z = 84 mG and for N = 6.8 × 10 4 atoms.
The above results motivate, and partially justify, a radical single-mode approximation (SMA) that assumes that the field operator can be written asψ α (r) ≈ Ψ 0 (r)â α , such that the espatial part can be integrated out and we can deal with the spin part in full. That is, the many-body aspects of the spin part can be fully taken into account. This will lead, as is the purpose of this paper, to show quantifiable aspects of intrinsic decoherence, as explained above, and the recurrence or revivals of the predictable oscillations for initial coherent states. The SMA approximation leads to a seemingly simple Hamiltonian for an Within SMA, this Hamiltonian differs from that of Eq. (1) by a term proportional to the number operatorN , which commutes withĤ. The first and second terms represent the linear and quadratic Zeeman contributions, p ∼pB, q ∼qB 2 , while the third one is the spin-mixing interaction, η ∼ (c 2 − c 0 ) such that η > 0 is polar and η < 0 ferromagnetic [32].
As described in the Introduction our interest here is the ferromagnetic case. From now on, our goal is the study of HamiltonianĤ and its ensuing dynamics and, therefore, we shall vary all three parameters independently in order to elucidate their role in the dynamics.
These parameters can certainly be tuned experimentally. The above operators inĤ are given in terms of the creation and annihilation operators of bosonic atoms in the z-direction spin states +1, 0 and −1, in obvious notation, Certainlyˆ f = (f x ,f y ,f z ) is the one-body spin or magnetization vector operator, witĥ andf 2 =f 2 x +f 2 y +f 2 z , and whose expectation value is the magnetization or spin texture. 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 (3) 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 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, M , n 0 |f 2 |M, n 0 = 2 (n −1 + 1)(n −1 + 1)n 0 (n 0 − 1)δ M M δ n 0 (n 0 −2) where n +1 = (N + M − n 0 )/2 and n −1 = (N − M − n 0 )/2. We note that the matrix elements in Eqs. 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 [46] and F = 2 [44]. The presence of the quadratic term, q = 0, breaks the axisymmetry [40] and this system becomes an excellent one to study Quantum Phase Transition [47], quench-dynamical behaviors [48], Quantum Kibble-Zurek Mechanism [49], spin fragmentation [50] 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. (7)- (9). Hence, reiterating, the presence of the atomic interactions is mediated by the presence of the Zeeman quadratic term.

ENT STATES
With the full quantum diagonalization described above we can calculate the time evolution of arbitrary initial states |Ψ 0 . Here, we will concentrate in a family of initial coherent states, introduced below. 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, where the 1/N factor is introduced such that the trace of the reduced density is always unity. This yields in turn that the expectation values of spin operators are also bounded by one. 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 [51]. A very interesting and useful property of the one-body density matrix given above, in the representation of theâ α , α = −1, 0, +1 corresponding to the z−direction spin basis, is that the one-body reduced density matrix is diagonal in such a basis, While we have not proved this result, we have extensively verified it and we believe it follows from the commutation off z with the Hamiltonian. The fact that the true stationary density matrix ρ (s) jk is diagonal in this basis, will greatly facilitate the elucidation of the (quasi) stationary states reached in the time evolution of an initial coherent state.
As we can accurately calculate the matrix ρ kj (t) for any time using Eq.(13), we now turn our attention to the initial set of coherent states, [35] where ζ 1 = e −iϕ cos θ A very important property to take into account is that these coherent states are eigenstates off 2 , that is,f 2 |θ, ϕ = N (N − 1)|θ, ϕ . We point out that the distribution of energy eigenstates in an arbitrary coherent state involve, in general, several if not many blocks of different values of the quantum number M . 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 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 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 ρ jk (t), given by (13), can be first expressed as, where U QI (t, 0) = e −i(qQ+ηf 2 )t/ (18) 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. 4 we show sequences of decoherences and recurrences in the weak N |η|/q 1, strong q/N |η| 1 and crossover N |η|/q ∼ 1 regimes. 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.

IV. DECOHERENCE AND RECURRENCES IN THE STRONG AND WEAK IN-TERACTING REGIMES
To proceed with the finding of the weak and strong interacting extremes, we note that the density matrix, Eq. (17), can be written in two alternative forms, |θ e i((k−j)pt/ +ϕ) e iδ k,j±1 (δ j,2 −δ k,2 )qt/ (19) or whereV and e T stands for the time-ordered exponential. The first form, Eq. (19), is appropriate for the study of the weak interaction N |η|/q 1, while the second one, Eq. (20), for the strong one q/N |η| 1. Note that the exact limits η = 0 and q = 0 are explicitly recovered in Eqs. (19) and (20).
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.
[36] and [37], that for usual harmonic oscillator coherent states [52] show collapses whose envelope is of a gaussian shape followed by revivals, with characteristic times that depend exclusively on the state and the Hamiltonian parameters. Furthermore, the recurrences are observed to appear at periodic times. Based on this and on the numerical evidence, we proceed now to present the heuristic (semi) analytic forms of the density matrix in the two extremes.
A. Strong interaction regime q/N |η| 1 In this regime the leading terms depend mostly on the quadratic Zeeman coupling q and the one-body reduced density matrix can be very precisely fitted by the following form, with ρ jk (0) given by, which was obtained using the lie algebra method [53] From Eq. (22, it can immediately be seen that the structure of periodic recurrences with decoherences in relatively shorter times than the former is given by the real exponential term exp N 2 sin 2 θ cos q(j−k) N 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 a gaussian function of the form e −t 2 /2τ 2 dec , from which we identify, Note that both times τ dec grow as does the number N of atoms that, as we discuss in the next section, appear as an "expected" feature of thermodynamic behavior in the limit of large number of atoms.
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, with a jj (θ) a term that it is at most of the order of 10 −9 , compared with the initial value ρ kj (0), for N < ∼ 10 3 . The decoherence time scales as τ termQ, remain essentially constant, the x and y spin components decohere to zero. Note from the figure that if q 1, but η finite, the recurrences can be separated by very long times that increase as N grows, thus yielding a true stationary state in the thermodynamic limit. To enquire into the nature of the stationary state, as illustrated in Figs. 5 and 6 and in accord with Eqs. (22) and (26), we observe that the (quasi) stationary one-body density matrix becomes diagonal with their diagonal terms numerically very close to the initial ones.
As suggested by ETH, we compute the average energy E = θ, ϕ|Ĥ|θ, ϕ and find the closest eigenenergy E M,m M of the system to such an average. We have found in all studied cases for this regime that the corresponding stationary density matrix ρ jk equals the density matrix ρ (s) jk of the corresponding energy eigenstate |M, m M , see Eq. (14). Although we do not claim that the stationary state is a thermal, this result agrees with the consequences of the Eigenstate Thermalization Hypothesis. That is, all the one-body properties in the (quasi) stationary states are the same as if the system were in a true, exact, stationary energy eigenstate with the same eigenvalue as the average energy of the time-evolving system. As depicted in Fig. 5, the recurrences become more separated as N increases, thus indicating that for a very large system, the quasi stationary state cannot be distinguished from a true eigenstate, as far as measurable few-body properties are concerned. This is also along the explanation of thermalization in isolated many-body systems in statistical physics [18]. ρ 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(2fjkηt)−1) for all j, k .
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. 8 and 9. The diagonal terms g jj (θ) are all equal, shown in Fig. 9. The initial condition is given again by the matrix in Eq. (23).
Figs. 7 and 10 show the density matrix ρ jk (t) and the expectation value of the operator f x for a single example in the limit N |η|/q 1 (η = −10 −5 and q = 7, for N = 300).  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) ≈ e g jk (θ)N (cos(2fjkηt)−1) within the density matrix, see Eq. (27). 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, e g jk (θ)N (cos(2fjkηt)−1) ≈ e −g jk (θ)N (2fjkηt) 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 that is very different to the strong interaction case, see Eqs. (24) and (25). This is further discussed in the last section.
With regard to the (quasi) stationary states reached between consecutive recurrences or revivals, as shown in Figs. 7 and 10 and verified in expressions given by Eqs. (27)- (29), the off-diagonal terms of the reduced density matrix again vanish for time intervals during the (quasi) stationary state, with the diagonal terms becoming equal to the initial state. As in the strong interaction regime, since the energy eigenstates one-body density matrices, see Eq. (14), are also diagonal in the same spin basis, this suggests the same test of comparing the density matrices in the stationary regime of the coherent states with that of a true stationary eigenstate with the same average energy. Interestingly, contrary to the strong interacting regime, in this case those do not agree. That is, the stationary state between recurrences does not seem to be a "true" stationary state, in the sense of ETH. This suggests, along the mentioned N dependence of the decoherence and recurrence times, that the weak regime is not typical of thermodynamic states. Although the full analysis of these stationary states are out of the scope of the present article, this indicates that the distribution of energy eigenstates in coherent states may show a more complex structure in the weak than in the strong limit.

V. DISCUSSION AND FINAL REMARKS
The agreement between the proposed heuristic expressions for the time evolution of the rec ∼ N . It is of interest to mention that the nonlinear single-mode approximation, as described in Refs. [36] and [37], corresponds to the weak-interaction regime. We find relevant to recall that the validity of thermodynamics is achieved for macroscopic bodies N 1. In this limit, averages of extensive quantities scale as N while their fluctuations as √ N such that the ratio of deviations from equilibrium values scale as 1/ √ N , hence tending to zero as N increases [18]. Moreover, the larger the system, the longer it takes to equilibrate in a stationary state and the recurrences to initial states also become apart for a longer time. In this sense, it appears that the strong interacting regime fulfills this limit. That is, τ rec ∼ N appear as reasonable results. It is the more appealing the fact that the ratio τ dec /τ rec ∼ 1/ √ N , the mentioned typical condition for thermodynamic stability of a macroscopic system [18]. Thus, 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. As we have also seen, this strong difference appears to be present in the structure of the stationary states attained between consecutive recurrences, having started in coherent states. As we have verified, those stationary states are indistinguishably from true stationary energy eigenstates. The result that this property does not hold for a weakly interacting SBEC, appears to be in agreement with the N dependence of the decoherence and recurrence times mentioned above. The full elucidation of these differences certainly deserves a separate and detailed study. In the crossover regime N |η| ∼ q both the interaction η and non-linear Zeeman q contributions compete, their effect is intertwined and even though the system indeed shows decoherences and recurrences, the latter no longer occur at prescribed times depending on either √ N or N , see Fig. 4 (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. A natural extension of the present study is the analysis 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, to find out if the predicted decoherence and recurrence time hold and it would also serve to further enquire into the structure and properties of the stationary state between recurrences.