Quantum Phases of Time Order in Many-Body Ground States

Understanding phases of matter is of both fundamental and practical importance. Prior to the widespread appreciation and acceptance of topological order, the paradigm of spontaneous symmetry breaking, formulated along the Landau–Ginzburg–Wilson (LGW) dogma, is central to understanding phases associated with order parameters of distinct symmetries and transitions between phases. This work proposes to identify ground-state phases of the quantum many-body system in terms of time order, which is operationally defined by the appearance of the non-trivial temporal structure in the two-time auto-correlation function of a symmetry operator (order parameter) while the system approaches thermodynamic limit. As a special case, the (symmetry protected) time crystalline order phase detects continuous time crystal (CTC). We originally discover the physical meaning of CTC’s characteristic period and amplitude. Time order phase diagrams for spin-1 atomic Bose–Einstein condensate (BEC) and quantum Rabi model are fully worked out. In addition to time-crystalline order, the intriguing phase of time-functional order is discussed in two non-Hermitian interacting spin models.


INTRODUCTION
A consistent theme for studying the many-body system, particularly in condensed matter physics, concerns the classification of phases and their associated phase transitions [20,52,68]. In the celebrated Landau-Ginzburg-Wilson (LGW) paradigm [35,70], spontaneous symmetry breaking plays a central role with order parameters characterizing different phases of matter possessing respective broken symmetries. Other schemes for classifying phases as well as their associated transitions are, however, beyond the Landau-Ginzburg-Wilson paradigm, which are by now well accepted since first established decades ago [53,63,64]. For example, topological order, which classifies the gapped quantum many-body system, constitutes a topical research direction [63,64,66,67]. Our current understanding categorizes gapped systems into gapped liquid phases [74] and gapped non-liquid phases, with the former broadly including phases of topological order [63,64], symmetry-enriched topological order [9,12,25,65], and symmetry-protected trivial order [10,11,23], while the recently discussed fracton phases [55,56,60] belong to the latter of gapped non-liquid phases.
Temporal properties of phases are also worthy of investigations as exemplified by many recent studies [41,50,69]. For instance, time crystal (TC) or perpetual temporal dependence in a manybody ground state that breaks spontaneously time translation symmetry (TTS) constitutes an exciting new phenomenon. First proposed by Wilczek [69] for quantum systems and followed by Shapere and Wilczek [54] for classical systems in 2012, TC in their original sense is unfortunately ruled out by Bruno's no-go theorem the following year [3,42]. Watanabe and Oshikawa (WO) reformulate the idea of quantum TC, and present a refined no-go theorem for many-body systems without too long-range interactions [62]. Continued efforts are directed at searching for continuous time crystal in open systems [4,5,30] and classical driven-diffusive systems [28]. Most recent efforts on this topic are directed toward non-equilibrium discrete/ Floquet TC breaking discrete TTS [16,17,31,51,61,72], particularly in systems with disorder that facilitate many-body localizations [61,72], in addition to clean systems [19,27,39,49]. Ongoing studies are further extended to open systems with Floquet driving in the presence of dissipation [14,21,22,37,46], with experimental investigations reported for a variety of systems [2,13,43,47,48,57,75]. A recent study addresses TC and its associated physics along the imaginary time axis [6].
We introduce time order in this work as the essential element for a new perspective to identify and categorize quantum manybody phases, based on different ground-state temporal patterns. Each quantum many-body HamiltonianĤ comes with its evolution or time translation operator e −iĤt . When continuous time translation symmetry is broken for operator e −iĤt , akin to the breaking of continuous spatial translation symmetry for operator e −i k· r , time crystals arise in direct analogy to spatial crystals [69]. The message we hope to convey here in this study is rooted on the dual betweenĤ and e −iĤt , which we argue quite generally establishes a solid foundation for time order and provides further information concerning ground-state quantum phases based on time domain properties. Different quantum many-body states with the same temporal patterns are classified into the same time order phases, of which continuous TC (CTC), a ground state with periodic time dependence breaking continuous TTS as originally proposed in Refs. [54,69], belongs to one of them.
We will adopt the WO definition of CTC based on two-time auto-correlation function of an operator. It was first outlined in the now famous no-go theorem work [62], and it establishes a general and rigorous subtype of CTC. Recently, Kozin and Kyriienko claim to have realized such a genuine ground-state CTC in a multi-spin model with long-range interaction [33], buttressing much confidence to the search for exotic CTCs. The operational definition that we introduced for time order encompasses WO CTC as one type of time order phases. We will also explore and elaborate a variety of possible exotic phases.

Time Order
We argue that ground-state temporal properties of a quantum many-body system can be used to characterize or classify its phases. Hence, the concept of time order can be introduced analogous to an order parameter by bestowing it in the nontrivial temporal dependence. To exemplify the essence of the associated physics, we shall present an operational definition for time order and accordingly work out the exhaustive list of all allowed phases. According to the WO proposal [62], a witness to CTC is the following two-time (or unequal time) auto-correlation function (with respect to the ground state): for operatorΦ(t) ≡ V d D xφ( x, t) defined as an integrated order parameter (over D-spatial-dimension), or analogously the volume averaged one, withφ( x, t) the corresponding local order parameter density operatorφ ≡Φ/V. If f(t) is time periodic in the thermodynamic limit, the system is in a state of CTC. This can be reformulated into an explicit operational protocol by introducing a twisted vector. For a quantum many-body system with energy eigen-state |ψ i 〉, if there exists a coarse-grained Hermitian order parameterφ, ϕ|ψ i 〉 is called the eigen-state twisted vector; more generally, if ϕ is non-Hermitian,φ|ψ i 〉 (orφ † |ψ i 〉) will be called the right (or left) eigen-state twisted vector. The orthonormal set of eigen-wavefunctions |ψ i 〉 (i 0, 1, 2, / ) for a system described by HamiltonianĤ is arranged in increasing eigen-energies ϵ i with i = 0 denoting the ground state. When the coarse-grained order parameterφ is Hermitian, the ground-state twisted vector |v〉 can be expanded |v〉 ≡φ(0)|ψ 0 〉 ∞ i 0 a i |ψ i 〉 into the eigen-basis. With the help of the Schrödinger equation iz|ψ(t)〉/zt Ĥ |ψ(t)〉 (Z = 1 assumed throughout) for the system wave function |ψ(t)〉, we obtain the following equation: where η j ≡|a j | 2 denotes weights of the ground-state twisted vector, η 0 the corresponding ground-state weight, and η j (with j > 0) the excited-state weight.
When the coarse-grained order parameterφ is non-Hermitian, we use |v (l) 〉 and |v (r) 〉 to denote, respectively, the left and right ground-state twisted vectors and expand them analogously in the eigen-basis to arrive at |v (l) 〉 ≡φ In this case, we find the following equation: with η j ≡ b p j a j weights of the ground-state twisted vector instead. Similarly, η 0 and η j (j > 0) denote, respectively, ground-and excited-state weights.
Given an order parameterφ, quite generally f(t) is a sum of many harmonic functions with amplitudes η j and characteristic frequencies ω j ≡ϵ j − ϵ 0 . Non-trivial time dependence of the twotime auto-correlation function is thus imbedded in the energy spectra of H as well as in the weights of the ground-state twisted vector. For CTC order to exist, one of the excited-state weights must be non-vanishing in the thermodynamic limit, or in rare cases, f(t) can include harmonic terms of commensurate frequencies.
If f(t) is a constant, the time dependence will be trivial. However, a subtlety appears when f(t) is vanishingly small with respect to the system size. Since what we are after is the system's explicit temporal behavior or time dependence, it is easily washed out to f(t) = 0 by a vanishing norm of the twisted vector. Such a difficulty can be mitigated by multiplying system volume V, that is, using the twisted vector |v〉 → V|v〉, to check if the correlation for the bulk order parameter F(t) ≡ V 2 f(t) exhibits temporal dependence, or vanishes as follows: When f(t) = 0 but F(t) remains a periodic function, the system can still be considered a CTC. Such a remedy surprisingly captures the essence of generalized CTC of Ref. [40]. The analysis presented above can be directly extended to excited states [59]. It is also straightforwardly applicable to non-Hermitian systems, as long as a plausible "ground state" can be identified, for example, by requiring its eigen-energy to possess the largest imaginary part or the smallest norm. Denoting the imaginary part of energy eigen-value E i as Im(E i ), a prefactor ∝ e Im(Ei)t then arises in the autocorrelation function, leading to unusual time functional order in the classification of time order.
Therefore, quantum many-body phases can be classified according to time order. The two-time auto-correlation function-based complete operational procedure for classifying time order thus extends the definition of WO CTC as provided in Ref. [62]. Our central results can be simply stated as follows: if f(t) exhibits non-trivial time dependence, then time order exists. If f(t) = 0 but F(t) displays non-trivial time dependence instead, then generalized time order exists.
More specifically, if f(t) = const. is non-zero, the system exhibits time trivial order. The same applies when f(t) = 0 and F(t) = const. For all other situations, non-trivial time order prevails. A complete classification for all time order groundstate phases is shown in Table 1, according to the temporal behaviors of their auto-correlation functions f(t) or F(t). As discussed in Section 4, the above discussion and classification on time order can be extended to finite temperature systems as well.
The operational procedure outlined previously presents a straightforward approach for detecting time order, albeit with reference to an order parameter operator. Hence, more appropriately, this approach should be called order parameter assisted time order or symmetry-based (or -protected) time order to emphasize its reference to symmetry order parameter of a quantum many-body system. The twisted vector facilitates easy calculations to distinguish between different time order phases from time trivial ones, as we illustrate in the following text in terms of a few concrete examples. It is reasonable to expect that transitions between different time order phases can occur, reminiscent of phase transitions in the LGW spontaneous symmetry breaking paradigm.

Time Order Phase in a Spin-1 Atomic Condensate
A spin-1 atomic Bose-Einstein condensate (BEC) under single spatial mode approximation (SMA) [36,44,73] is described by the following Hamiltonian: whereâ mF (m F 0, ± 1) (â † m F ) denotes the annihilation (creation) operator for atom in the ground-state Zeeman manifold |F = 1, m F 〉 with corresponding number operator N m F â † m Fâ m F . The total atom numberN N 1 +N 0 +N −1 is conserved. p and q are linear and quadratic Zeeman shifts that can be tuned independently [38], while c 2 describes the strength of spin exchange interaction.
The validity of this model is well established based on extensive theoretical [8,24,71,76] and experimental [1,7,38,45] studies of spinor BEC over the years. The fractional population in spin states |1, 1〉 and |1, − 1〉,n sum ≡ N sum /N, with N sum N 1 +N −1 N − N 0 , is often chosen as an order parameter [1,15,34,71] with N assuming the role of system size. The ground state twisted vector then becomes |v〉 ≡n sum |ψ 0 〉, and TABLE 1 | Classification of the ground-state phases for a quantum many-body system.

Phase
Property of two-time auto-correlator Time trivial order

Time order
Time crystalline order f(t) is periodic and non-vanishing Time quasi-crystalline order f(t) is quasiperiodic with beats from two incommensurate frequencies Time functional order f(t) is aperiodic Generalized time crystalline order f(t) = 0, F(t) is periodic and non-vanishing Generalized time quasi-crystalline order f(t) = 0, F(t) is quasiperiodic with beats from two incommensurate frequencies Generalized time functional order We will concentrate on the zero magnetization F z = 0 subspace and employ exact diagonalization (ED) to calculate eigen-states. p = 0 is assumed since F z is conserved. Figure 1 illustrates the system's complete time order phase diagram. For ferromagnetic interaction c 2 < 0 as with 87 Rb atoms, the critical quadratic Zeeman shift q/|c 2 | = 2 splits the whole region into the time trivial order (TT) phase for smaller q that observes TTS and the generalized time crystalline (gTC) order phase for q/|c 2 | > 2, where TTS is spontaneously broken. The latter (gTC phase) is found to coincide with the polar phase [71]. Limited by available computation resources, the system sizes we explored with ED remain moderate which prevent us from mapping out the finer details in the immediate neighborhood of q = 2|c 2 |. Further elaboration of time order properties in this region is therefore needed. On the other hand, for antiferromagnetic interaction c 2 > 0 with 23 Na atoms, we find q = 0 separates TT phase from gTC order. We note here that q = 2|c 2 | is the second-order quantum phase transition (QPT) critical point between the polar phase and the broken-axisymmetry phase of the ferromagnetic spin-1 BEC, while q = 0 corresponds to the first-order QPT critical point for antiferromagnetic interaction.
More detailed discussions including the dependence of time order phases on system size, possible approaches to detect them, and extension to thermal state phases can be found in Section 4.

Time Order Phase Diagram for Quantum Rabi Model
As a second example, we consider time order phases of the quantum Rabi model described by the Hamiltonian as follows: whereσ x,z is the Pauli matrix of a two-level system (transition frequency Ω),â(â † ) is the annihilation (creation) operator for a single bosonic field mode (of frequency ω 0 ), and λ is their coupling strength.
It is known that the aforementioned model exhibits a QPT to a superradiant state, despite its simplicity [29]. The transition occurs at the critical point g c ≡ 1, with the dimensionless parameter g ≡ 2λ/ ω 0 Ω √ . The equivalent thermodynamic limit is approached by taking Ω/ω 0 → ∞. Though the system only has finite components, the QPT herein is well established. According to the studies in Ref. [29], an almost exact effective low-energy Hamiltonian for the normal phase (g < 1) is given by the following whose low-energy eigen-states are |ϕ m For the superradiant phase (g > 1), the effective low energy Hamiltonian becomeŝ whose eigen-states are given by More details can be found in the supplementary material of Ref. [29].
For this model, the scaled average cavity photon numbern c ω 0â †â /Ω is a suitable order parameter with Ω/ω 0 assuming the role of system size. The corresponding bulk order parameter then becomesN c â †â or the average cavity photon number, and For g < 1, we find respectively, where η 0 = sinh 4 (r np ) and η 2 = cosh 2 (r np ) sinh 2 (r np ). For g > 1, we obtain the following equation: The time order phase diagram is shown in Figure 2. When g < 1, the system ground state corresponds to a generalized time crystalline order phase, while the system exhibits time trivial order when g > 1. Despites such a simple model composed of a two-level system and a bosonic field mode, the ground state of the quantum Rabi model displays an intriguing temporal phase structure accompanied by a finite-component quantum phase transition.

Non-Hermitian Many-Body Interaction Model
Finally, we consider two effective models with many-body spin-spin interaction and non-Hermitian effects. The first is described by the Hamiltonian: with two σ y operators at sites i and j in a string of otherwise σ x Nbody spin interaction. 1/[N(N − 1)] is the normalization factor, λ is the spin interaction strength, and γ represents an effective dissipation rate. λ > 0 and γ ≥ 0 are both real numbers. We observe that the Greenberger-Horne-Zeilinger (GHZ) states correspond to two non-degenerate system eigen-states with eigen-energies ± (λ + iγ)/2. The spectra of this model system are bounded inside the circle of radius λ 2 + γ 2 /2 in the complex plane. The eigen-state whose eigen-value has the largest imaginary part is taken as the ground state, or |GS〉 = |G + 〉 with eigen-energy ϵ 0 = (λ + iγ)/2. The highest excited state is |G − 〉, whose corresponding eigen-energy is ϵ 2 N −1 −(λ + iγ)/2. An appropriate order parameter operator in this case becomes the average magnetizationm N i 1 σ z i /N. The twisted vector becomes |v〉 m|GS〉 |G − 〉, and the auto-correlator can be easily worked out to be f(t) lim N→∞ 〈m(t)m(0)〉 e −iλt e γt . When γ = 0, the system ground state exists time-crystalline order phase and corresponds to a continuous time crystal [33]. When γ ≠ 0, the system exhibits time functional order, with an exploding f(t) as time evolves.
Therefore, we directly obtain the following equation: When λ ≠ 0 and γ ≠ 0, the system exists in time functional order phase and again results from the non-Hermitian Hamiltonian. When λ ≠ 0 but γ = 0, the auto-correlation function reduces to as for a genuine time crystal of the WO type exhibiting time crystalline order. When λ = 0 and 0 < |γ| ≤ 1, we find The system ground state again exhibits time-crystalline order. When λ = 0 and |γ| > 1, we obtain by choosing ϵ 0 −N + 2 + 2i γ 2 − 1 as the ground state eigenenergy from the two eigen-values −N + 2 ± 2i γ 2 − 1 . The system ground state now exhibits time functional order phase, with a decaying f(t) as time evolves. When λ = γ = 0, the ground state reduces to the time trivial order phase. The above two non-Hermitian models represent direct generalizations of the Hermitian system as considered in Refs. [18,33]. While slightly more complicated, they remain sufficiently simple for compact analytical treatment, thus helping to reveal interesting and clear physical meanings of the underlying time order.

Some Remarks About Continuous Time
Crystal exhibits no temporal dependence; hence, it belongs to time trivial order according to our classification scheme. At first sight, this seems to sweep many important models of condensed matter physics into the same boring class of the time trivial order phase. However, it remains to explore, for instance, many-body systems with more than two-body (or kbody) interactions, or non-Hermitian systems, which might support the existence of CTC. Inspired by the recent results on CTC [33], we believe more time crystalline phases will be uncovered and further understanding will be gained in the future.
The two-time auto-correlation function f(t) measures the CTC phase as in earlier studies [33,62], while both autocorrelation functions f(t) and F(t) together define different time-order phases we propose in this work. The absence of the local temporal behavior f(t) = 0 does not imply the absence of any temporal property in the bulk, when F(t) could have various temporal behaviors. Based on this, our operational definitions of time order are developed. We also extend the scope of CTC to include both TC order and gTC order phases. This distinction between f(t) and F(t) gives more insights into quantum many-body phases.
As emphasized earlier, continuous time crystal results from spontaneously breaking continuous time translation symmetry. Due to the genuine time periodicity contained in CTC, it might be possible to explore and design new types of clocks based on macroscopic many-body systems, as the time period is directly related to energy spectra, and whose physical meaning is clearly the same as for atomic clock states. Furthermore, they are not affected by finite size effect in contrast to periodicity in DTC.

DISCUSSION AND CONCLUSION
While ground-state phases of a quantum many-body system are mostly classified with their Hamiltonian based on the following two paradigms: LGW symmetry breaking order parameter or topological order, this work proposes to study phases from time dimension using time order or more specifically with the proposed symmetry-based time order. Compared to the recent progress and understanding gained for topological order [66,67], one could try to develop a framework for entanglement-based time order instead of the symmetry-based time order we employ here in this study. Quantum entanglement in a many-body system is responsible for topological order, whose origin lies at the tensor product structure of the quantum many-body Hilbert space H tot ⊗ i H i with H i the finite-dimensional Hilbert space for site-i. An entanglement-based time order therefore calls for a combined investigation to exploit quantum entanglement and temporal properties of a quantum many-body system.
Through time order, one focuses on temporal structure of the evolution operator e −iĤt . The symmetry-based time order therefore unifies the LGW paradigm with the concept of time order, while an entanglement-based time order could amalgamate the topological order paradigm (or entanglement beyond that) with time order. For this to happen, a more basic definition for time order will be required, which will likely expand into further in-depth investigations.
In conclusion, understanding the phases of matter constitutes a cornerstone of contemporary physics. Capitalizing on the concept of CTC for the many-body ground state with perpetual time dependence, this study argues that information from time domain can be employed to classify the quantum phase as well, which provides a new perspective toward the understanding of ground-state time dependence, significantly beyond existing studies on CTC. We introduce time order, provide its operational definition in terms of two-time autocorrelation function of an appropriate symmetry order operator, bestow physical meaning to characteristic frequencies and amplitudes of the correlation function, and present a complete classification of time order phases. Time order phase diagrams for a spin-1 BEC system and the quantum Rabi model are fully worked out. Interesting time order phases in non-Hermitian spin models with multi-body interaction are presented. In addition to the time crystalline order which already attracts broad attention from its studies in terms of CTC, other phases we identify, for example, time quasi-crystalline order and time functional order, represent exciting new possibilities.

METHODS AND CALCULATION DETAILS
Here, in this section, we provide more supporting material for our main results and related details for the aforementioned presentation. It is organized as follows: in Section 4.1, we extend the discussion of time order to finite temperature; in Section 4.2, we present calculation details related to the spin-1 atomic BEC example considered and give intriguing results for finite temperature scenario in spin-1 BEC; in Section 4.3,as a more straightforward approach to understand numerical results, we present a variational approach for treating the polar ground state of a spin-1 BEC. Finally, we give the details about ground state and eigen-energy calculation in the non-Hermitian quantum many-body model with multi-body interaction in Section 4.4.

Time Order at Finite Temperature
At finite temperature T, excited states will be populated, which can be taken into account with the Gibbs ensembleρ ≡ e −βĤ /Z, where Z ≡ Tr e −βĤ denotes the partition function and β ≡ 1/T the inverse temperature. We then find where |v k 〉 is the eigen-state twisted vector for |ψ k 〉 and c jk is its associated weight. Analogously, for the non-Hermitian case, we find where |v (l) k 〉 and |v (r) k 〉 are the left and right twisted vectors for eigen-state |ψ k 〉, respectively, and c jk is the corresponding weight.
It is easily noted that f(t) at finite temperature contains contributions from all eigen-states of the quantum many-body systemĤ, with a temperature-dependent weight factor for different energy levels, but f(t) remains to include contributions from different periodic functions. Hence, the quantum phase classification task essentially remains the same (including its possible reference to F(t)) as is shown in Table 1 for the ground state. At finite temperature, due to thermal excitations to the ground state, the temporal behavior will be more complex, thus opening up for more interesting possibilities, for example, to control time order phases and to study crossover or driven phase transitions between different time order phases.

Time Order in a Spin-1 Atomic BEC
For typical interaction parameters of a spin-1 BEC (e.g., of ground state 87 Rb or 23 Na atoms) in a tight trap, spin domain formation is energetically suppressed when the atom number is not too large as spin-dependent interaction strength is much weaker than spin-independent interaction [26,32,36,58]. This facilitates a single-spatial-mode approximation (SMA) by assuming all spin states share the same spatial wave function ϕ(r), which effectively decouples the spatial degrees of freedom from the spin and results in the Hamiltonian [36,44] in Eq. 6 for the model many-body system, whereâ mF (m F 0, ± 1) is the annihilation operator of the ground manifold state |F = 1, m F 〉 with corresponding number operatorN m F â † m Fâ m F . p and q are linear and quadratic Zeeman shifts which could be tuned independently in experiments [38], while c 2 denotes the spin exchange interaction strength. Unless otherwise noted, we will take |c 2 | = 1 as unit of energy in this work. The total particle number operatorN N 1 +N 0 +N −1 and the longitudinal magnetization operatorF z N 1 −N −1 are both conserved. Thus, linear Zeeman shift can be set to p = 0 effectively.
As discussed in the main text, a suitable order parameter for this model system isn sum ≡N sum /N (N sum N 1 +N −1 N −N 0 ), which measures the fractional atomic population in the states |1, 1〉 and |1, − 1〉, and N assumes the role of system size. Following our formulation and denoting the system energy eigen-state by |ψ i 〉 (i = 0, 1, 2, /) with increasing eigen-energy ϵ i , the ground-state twisted vector becomes |v〉 ≡n sum |ψ 0 〉 ∞ i 0 a i |ψ i 〉, with a i = 〈ψ i |v〉 its expansion coefficient on the eigen-state |ψ i 〉. We find where b j ≡|a j | 2 is the weight of the ground-state twisted vector, b ≡ ∞ j 0 b j the total weight, and where A i = N〈ψ i |v〉, B j ≡|A j | 2 is the weight of the enlarged ground state twisted vector, and B ≡ ∞ j 0 B j the total weight. Our study below is for the zero magnetization F z = 0 subspace and employs exact diagonalization (ED) to calculate eigen-states as well as eigen-energies. The overall time order phase diagram for spin-1 BEC is shown in Figure 1. For ferromagnetic interaction c 2 < 0, the critical quadratic Zeeman energy q/|c 2 | = 2 splits the whole region into the time trivial order (TT) phase for smaller q that observes TTS, and the generalized time crystalline (gTC) order phase for q/|c 2 | > 2 where TTS is spontaneously broken. The latter (gTC phase) is found to coincide with the ground-state polar phase. The available computation resource limits the calculation to a finite system size, which prevents us from mapping out the exact details in the immediate neighborhood of q = 2, where further elaboration is needed for its time order properties. On the other hand, for antiferromagnetic interactions, we find q = 0 separates TT phase from gTC order.
In Figure 3, the weights for the ground state as well as for the low-lying excited states are shown as functions of q for a typical system size of N = 10 000. Only the ground-state weight b 0 is nonvanishing in the q < 2 (q < 0) region for ferromagnetic (antiferromagnetic) interactions, but total weight b is zero in the q > 2 (q > 0) region for ferromagnetic (antiferromagnetic) interaction, which prompts us to examine further the enlarged weights B i corresponding to the bulk order parameter. For ground and the first excited states, the volume enlarged weights B 0,1 are found to be non-vanishing, although both decrease as q increases and grow with N as q approaches q = 2 (q = 0) for ferromagnetic (antiferromagnetic) interaction. However, as mentioned above, limited to a system size of N = 10 000 by computation resource in the ED calculation, we cannot exactly map out the behavior near q = 2 (q = 0) for ferromagnetic (antiferromagnetic) interaction. This consequently leaves empty for q in region [2.0, 2.02] ([0, 0.01]) for ferromagnetic (antiferromagnetic) interaction. The dependence on system size N is clearly revealed by Figure 4, with the enlarged weights in the gTC regime attaining fixed values as the system approaches thermodynamic limit (N → ∞). In regions away from q = 2 (q = 0) for ferromagnetic (antiferromagnetic) interaction, ED numerics can always approach thermodynamic limit, except for the immediate neighborhood near q = 2 (q = 0), where we infer with confidence the tendencies to divergence of the weights B 0,1 as q approaches q = 2 (q = 0).
The time evolution of two-time auto-correlation function F(t) is plotted in Figures 5A,C for ferromagnetic and for antiferromagnetic interactions, while Figures 5B,D display energy gaps between ground and the first excited states as a function of q for ferromagnetic and antiferromagnetic interactions, respectively, at a system size of N = 5000. The behavior of F(t) is quantitatively consistent with that of the weights B i (q) (i = 0, 1) shown in Figure 3 and the energy gap ϵ 1 − ϵ 0 shown in Figures 5B,D. At finite temperature, excited states come into play by also contributing to the correlation function. We find the gTC order hosted in the polar phase persists for both ferromagnetic and antiferromagnetic interactions. The corresponding time evolution and Fourier transform of F(t) are shown in Figure 6, calculated for N = 500 at a temperature of β ≡ 1/T = 1. The Fourier transform is performed for Re(F) over t = [0, 1000] with the zero frequency (DC) component subtracted or for Im(F). The upper (lower) panel corresponds to ferromagnetic (antiferromagnetic) interaction at q = 3 (q = 2). For ferromagnetic interaction, two distinct frequency components are clearly identified for q = 3, associated with the two different energy level gaps. The beautiful beat pattern for  F(t) would appear, while we only show the short time behavior in Figure 6 (a). Thus, the gTC phase remains at a finite temperature. Moreover, we also find a generalized time quasi-crystalline order phase assuming the two frequencies are incommensurate, by finetuning their corresponding energy gaps such that the relation Δ 1 /Δ 2 = m 1 /m 2 with m 1 and m 2 being co-primes is not satisfied. The gTC  phase at finite temperature here is robust which is in contrast to the melting behavior of continuous time crystal (CTC) shown in Ref. [33]. Finally, we hope to address the critical question about how could this time order, sort of a perpetual time dependence, can be observed. We note the bulk two-time auto-correlation function introduced F(t) lim N→∞ 〈N sum (t)N sum (0)〉 denotes nothing but the ground-state (averaged) conditional outcome of measuring N sum (t) at t after starting with N sum (0) initially. The dynamics of F(t) follows that of N sum (t) as in the quantum regression theorem. Given the system is well controlled, highly reproducible, one can simply detect F(t) by measuring N sum (t), although for each measurement at an instant t, a condensate is destroyed, and a follow-up one will have to be prepared as closely as possible in every respects (through selection and post-selection) and be measured at a different t′ > t. Thus, a plausible way to detect the ground-state time dependence will require reconstructing the time dependence of F(t)/N sum (0). As long as the oscillation amplitude is more than a few percent, it will be easily observable with not too much difficulty, although such a reconstruction will still be difficult as N sum (0) can be rather small compared to N 0~N in the polar state. Alternatively, one can perhaps start from a twin-Fock state, that is, by preparing an initial state with N sum (0)~N.
In Figure 7A, we show the behavior of oscillation amplitude for F(t)/N sum (0). The time dependence of F(t)/N sum (0) at q = 2.5 for ferromagnetic interaction is shown in Figure 7B

A Variational Polar State for Ferromagnetic Spin-1 BEC
One might naively expect that nothing particularly interesting could happen in the polar phase of a ferromagnetic spin-1 BEC, where essentially all atoms reside in the single particle state |1, 0〉. Nevertheless, due to the competition between spin exchange interaction c 2 and quadratic Zeeman shift q, the ground state of our system differs from |N 1 = 0, N 0 = N, N −1 = 0〉, which can be affirmed based on a simple variational analytical calculation given in this section.
We use the number-state basis |N 1  is a (complex) variational parameter with r and ϕ as real parameters. From Eq. 6 and (assumed) p = 0, the ground-state energy follows from We see the extreme value (the minimum) of E is reached when cos(ϕ) = ±1, that is, for a real variational parameter a, which will be assumed from now on. This gives the following equation: which determines the locations for the extreme values.
and the corresponding extreme values are as follows: In the thermodynamic limit N → ∞, they reduce, respectively, to a ± 1 + q c2 ± 2c 2 2 +2c 2 q+q 2 √ c2 and E ± c 2 + q ± 2c 2 2 + 2c 2 q + q 2 . The left and right asymptotic value for the energy function E(a) is therefore For ferromagnetic interaction (c 2 < 0), E − assumes the minimum, which corresponds to the ground state |ψ 0 〉 Despite the vanishing order parameter n sum in the polar phase (here the gTC order phase from the time order perspective), the enlarged quantity N sum retains a finite value. Hence, the physics we present here clearly belongs to the realm of quantum effects, beyond the reach of mean-field theory.
Here, we choose the eigen-state from {|Ψ (+) 〉, |Ψ (−) 〉} as the ground state |Ψ 0 〉 of our generalized non-Hermitian system, for it deforms into the ground state of the Hermitian case when γ approaches zero. If ϵ is real, then ground-state energy ϵ 0 corresponds to the smaller one from ϵ (±) . However, ground-state energy ϵ 0 corresponds to the one with the larger imaginary part when ϵ is a complex number. Ground state |Ψ 0 〉 is obtained straightforwardly. For the GHZ state |G + 〉, we can know it is a non-degenerate excited state with energy ϵ + = − N, for

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
T-CG proposed and conducted the research, LY supervised the research, and T-CG, LY discussed the results and wrote the manuscript. ACKNOWLEDGMENTS T-CG thanks Qi Liu and Ming Xue for their helpful discussion about spinor Bose-Einstein condensate.