{'fraction_non_alphanumeric': 0.03481820370816277, 'fraction_numerical': 0.006790272092463279, 'mean_word_length': 4.776356050069541, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Imaginary time crystal of thermal quantum matter

Spontaneous symmetry breaking is responsible for rich quantum phenomena from crystalline structures to superconductivity. This concept was boldly extended to the breaking of time translation, opening an avenue to finding exotic phases of quantum matter with collective time modulation and correlation. Here we report that a thermally open quantum ensemble manifests in the dual space of imaginary time with crystalline ordering due to a bath-induced retarded interaction. Exact quantum Monte Carlo simulations are performed to show that this imaginary time crystal phase exhibits characteristic ground-state and thermal properties absent in conventional quantum manybody systems, especially the striking temperature-oscillating behavior of its physical observables.

Statistical field theory shows that a d-dimensional quantum thermal ensemble is equivalently described by the path integral representing the partition function in d + 1 Euclidean space where time is imaginary. In the field theory dual description, the thermal distribution weighted by the well-known temperature dependent factor is fully captured by integration over auxiliary imaginary time ("iTime") ranging from 0 to β = 1/k B T (k B the Boltzmann constant). In this manner, finitetemperature physics is mapped to the dynamics of Euclidean action. Motivated by such a profound relation between temperature T and iTime, Wilczek had speculated the novel possibility of spontaneous translational symmetry breaking in iTime in the end of his original paper, dubbed as "imaginary time crystal" (iTC) [1]. One may wonder what class of thermal ensembles can show crystalline order in iTime space, whether the "time" analogue of phonons exist and what novel effects should manifest in ground state and thermal response functions.
Here, we find a class of open quantum ensembles that shows the iTC phase. A key feature of our mechanism is a non-monotonic bath-induced retarded interaction, whose dependence is on the relative time in two body scattering. That satisfies the two conditions what previously might have seem conflicting, i.e., being time dependent but fully translational-invariant. As we will show in the following, integrating out certain engineered local thermal baths yields such a retarded interaction, which favors crystalline patterns in iTime with periods determined by the position of interacting potential minimum. Thermodynamic quantities of the iTC exhibit a striking behavior of oscillation in temperature, which provides an experimental observable effect for this novel phase.

II. MODEL AND METHOD.
We study hard-core bosonic models in one-dimensional (1D) lattice of L sites as well as a L×L square lattice with where J denotes the hopping amplitude and ij indicates a pair of adjacent sites. We focus on the half-filling case throughout this paper. On each site i the density operator of the hard-core boson n i additionally couples to a local bath, as shown in Fig. 1(a). We assume the total system (system+bath) is in thermodynamic equilibrium with temperature T . The partition function takes the form Z tot = Tr s Tr b e −βHtot (Tr s/b denotes tracing over the system/bath degree of freedom). We first integrate out the bath degrees of freedom, which yields a retarded density-density interaction term in iTime. Thus the open system can be described by the reduced density matrix where H s is the system Hamiltonian and Z s = Tr sρs . S R describes effective action of the translational-invariant onsite bath-induced retarded interaction: where n 0 = 1 2 is the average particle density, and the site-independent kernel function D(τ − τ ′ ) depends on the bath properties. For instance, a harmonic oscillator bath with frequency ω 0 gives rise to the kernel function We will show later that the iTC phase can exist for a broad class of kernel functions; here, we chose: where |τ | ≤ β, and α is the retarded interaction strength.
The model with the retarded interactions defined in Eq. (3) can be more than a toy model of merely academic interest. In general, the interaction induced by baths composed of free bosons in thermal equilibrium is usually attractive, while a fermionic bath-induced interaction is more complex as sometimes it can be oscillating between attraction and repulsion (e.g.the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction [21][22][23]). Now we need a RKKY-like interaction, but with oscillating decay in iTime instead of real space. Consider the Matsubara frequency and Σ b (iω m ) the self-energy of the bath. In general, an imaginary part in self-energy leads to an exponential decay of the quasi-particle in real time (finite lifetime), which corresponds to an oscillation in the imaginary time after analytic continuation. Motivated by Kozii and Fu's recent work about an effective non-Hermitian description for fermionic quasiparticles [24], In Appendix A, we propose a similar microscopic bath model that can induce the oscillating retarded interaction as shown in Eq. (3). Moreover, by analogy with interaction in conventional crystals (e.g. Van der Waals force), we should emphasize that the time-oscillating feature in the interaction is not a necessary condition for the iTC phase. In Appendix B, we have chosen a different form of retarded interaction and numerically shown that the phase can actually exist as long as the retarded interaction is nonlocal and non-monotonic in itime with at least one minimum, whose position defines the "time" lattice constant of iTC.
Path integral Quantum Monte Carlo (QMC) simulation is a stochastic numerical method to study the equilibrium properties of quantum many-body systems based on importance sampling of the world line configurations in Euclidean spacetime. Following the algorithm proposed in Ref. [31], we generalize the QMC simulation with worm-type updates [32,33] to study the open systems whose distribution functions [Eq. (1)] deviate from the Boltzmann distribution. The key point is to evaluate the integrals resulting from the retardation and include them into the QMC acceptance ratio during the updates of sampling. Notice that since the retarded interactions [Eq. (2)] only appear in the diagonal parts (under the Fock basis) of the effective actions, they do not cause extra minus sign problem. Thus our QMC simulations can be considered as numerically exact since our "system" Hamiltonian (of hard-core bosons) is also positive definite. Retarded interactions have also been implemented in other different QMC algorithms [34,35].
The bosonic models with retarded interactions can be solved numerical exactly by Quantum Monte Carlo simulations (see Methods). We first focus on the ground state of the total system (system+bath). In finite size scaling, the inverse temperature is chosen as β = L(2L) for 1D (square) lattice, and the ground state in the thermodynamic limit is approached in the limit L → ∞. In the 1D case without retardation, the ground state is a quasisuperfluid ground state characterized by algebraic decaying correlation functions in both space and iTime. As the retarded interaction is switched on, spatial-temporal configurations with periodically oscillating world lines are favored since the interacting energies are lowered in these patterns. The competition between the hopping and the retarded interactions will lead to novel quantum phases.
The crystalline order can be identified from the density correlation functions in space and iTime: S(r, τ ) = . We first focus on the onsite correlations with r = 0. As shown in Fig. 5(a), for small α, S(0, τ ) still decays algebraically with τ , but is accompanied by oscillations. While for sufficiently large α a long-range order is built up and a characteristic time period emerges, indicating a spontaneous breaking of translational symmetry in iTime, characterized by the order parameter m 2 τ = 1 β β 0 dτ S(0, τ ) cos 2πω c τ. A finite-size scaling of m τ is shown in Fig.5(c), from which we can find that for large α, m τ extrapolates to a finite value in the thermodynamic limit, indicating a crystalline order in iTime. The existence of a quantum phase transition can be further verified by the correlation length ξ τ along iTime (S(0, τ ) ∼ e −|τ |/ξτ for sufficiently large τ ), which can be derived from the structure factor (Fourier transformation of S(0, τ )) [25]. The normalized correlation length ξ τ /β as a function of α for different system size and β is plotted in Fig. 5(e), where we find a crossing point at α τ c = 0.23(1), indicating a scaling invariant critical point. In the iTC phase, ξ τ /β linearly scales with the system size, and extrapolates to a finite value in the thermodynamic limit (inset of Fig. 5(e)).
For sufficiently large α, it is interesting to notice that the crystalline pattern not only emerges in iTime, but also in real space, which can be identified by the equaltime correlation function S(r, 0) as well as its structure factor S(Q) = 1 L r e iQr S(r, 0). In the case of halffilling, the spatial crystalline pattern is a charge-densitywave (CDW) state with a periodicity twice that of the lattice, as shown in Fig. 5(b). A finite size scaling of the order parameter m s = S(Q = π) is shown in Fig. 5(d).
Since there is no direct nearest-neighboring repulsive interaction in our system Hamiltonian, the CDW order is formed by an effective repulsive interaction induced by the retardation: the retarded interaction (3) favors oscillating world lines in iTime, and thus requires the hardcore bosons being separated as far as possible to avoid blocking the oscillations of each other. Therefore at halffilling the tendency towards CDW order is expected in real space. Similar to the analysis above, we also plot the normalized correlation length ξ s /L in real space as function of α in Fig.5(f), where a continuous phase transition is found at α = 0.245 (5). Above this value, long-range CDW correlations emerge in real space.
To further investigate the properties of the iTC, we calculate the single particle Green's function: G(r, τ ) = a † i (τ )a i+r (0) . The equal-time part G(r, 0) is shown in Fig. 3(a), which decays algebraically with r (accompanied by an oscillation) even for large α, indicating a quasisuperfluid order in the iTC. This order can be further verified from the superfluid density, which can be calculated as ρ s = L β W 2 , where W is the winding number denoting the net times by which the world-lines wrap around the 1D lattice. As shown in the inset of Fig. 3(a), ρ s grows monotonically with α, indicating that the retarded interactions enhance the mobility of bosons. Now we turn to the temporal part of G(r, τ ), and focus on its k = 0 component: G(k = 0, τ ) = 1 L r G(r, τ ). Fig. 3(b) shows that G(k = 0, τ ) decays exponentially with τ , indicating a finite single-particle charge gap. This agrees with the results of the compressibility κ = β L ( N 2 − N 2 ) with N the total particle number. Fig. 3(c) show that κ decays to zero for large α, thus the iTC is incompressible.
In summary, the iTC phase induced by the retarded interaction (3) spontaneously breaks the translational symmetries in both real space and iTime. Besides these intriguing properties, this phase is unique in the sense that it is incompressible but possess quasi-long range superfluid order, so it is fundamentally different from the conventional bosonic ground states in any closed quantum many-body systems: e.g. the (quasi-)superfluidity, Mott insulator, supersolid, Bose glass or Bose metal. The nature of this phase can be qualitatively understood in the strong coupling limit α ≫ J, where the retarded interaction dominates in the partition function. For each site, the retarded interaction (3) favors configurations with alternating n(τ ) that oscillates between 0 and 1 with a frequency 2πω c , which is reminiscent of the soliton solutions in the BCS pairing model [26,27], similar states have recently been proved to be metastable in a wide class of closed quantum many-body systems [28]. Since a world line needs to be continuous, the temporal-spatial configurations minimizing the retarded interacting action (2) are highly degenerate. As shown in the Fig. 3(d), these configurations can be understood as a synchronous movement of the bosons in iTime, in the sense that different bosons hop simultaneously towards the same direction; These synchronous movements can reduce the possibility of the collisions of world lines of different hard-core bosons, and at the same time, minimize the retarded interacting energy. In this hugely degenerate manifold, those configurations with nonzero winding number contribute to the nonzero superfluid density (e.g. the 2nd figure in Fig.3(d)). The incompressible nature can also be understood in such a picture: an extra particle added into the half-filled system will inevitably block or impede the synchronous movement of other bosons due to its hard-core nature, and thus will increase the retarded interacting energy and give rise to a finite charge gap.
Symmetries and excitations. The iTC phase simultaneously breaks the continuous and discrete translational symmetry in temporal and spatial directions. The spacial and temporal spontaneous symmetry breakings lock with each other, leaving an unbroken subgroup of intertwined space-iTime translational symmetry. The unit cell of the iTC is spanned by the primitive basis vectors b 1/2 with a rhombic symmetry as shown in Fig. 3(e), as opposed to a "decoupled" direct product of vectors in τ and x directions with a rectangular symmetry. Due to the fact that the time continuum is locked with the discrete spatial lattice we started with, small translational fluctuations of the atomic positions along the primitive basis vectors do not generate gapless "phonons", but excitations with finite energy cost due to the discreteness of the spatial translational symmetry in our model. The space-iTime locking avoids the dangerous gapless modes in renormalization, so remarkably stabilizes the iTC phase.
IV. TEMPERATURE EFFECT OF ITIME CRYSTAL.
The most striking effect of iTC phase can be found at finite temperature. In conventional systems, the temperature dependence of a physical quantity could be either monotonic or non-monotonic with one or several peaks. In contrast, an iTC could exhibit exotic thermal behav-ior that its physical observable depends on temperature in an oscillating way. This intriguing property can be understood as follows: the crystal structure in iTime will perfectly fit at certain temperatures satisfying the condition that β is a multiple of the period of iTC; otherwise, for general β the mismatch between them will introduce defects or distortions, which impede the synchronous movement of the bosons, thus are detrimental to crystalline order in both iTime and space. As an example, in Fig. 3(f), we plot the CDW ordering parameter as a function of β, which shows an oscillating behavior with peaks at the positions of integer multiples of the iTC period. This abnormal thermal behavior can be considered as an experimental diagnostic of the iTC phase. Extension to square lattices. Up to now, we have focused on 1D in the spatial part, while a generalization of our method to 2D square lattices is straightforward. Similar to the 1D case, we find that the temporal crystalline order also emerges for sufficiently large retarded interactions (left inset of Fig. 4). For a fixed L and β there seems to be a jump of the iTC order parameters m 2 τ as a function of α (Fig. 4), indicating that it is likely a 1st order quantum phase transition. However, due to the finite size effect and the limited system size in our simulation, we cannot preclude other possibilities. Another important difference between 1D and 2D is that in the latter case, there exists true superfluidity long-range order characterized by the condensate fraction ρ c = 1 L 2 r G(r, 0). In Fig. 4 (b), we plot ρ c as a function of α, where we also find a discontinuity occurs at α c = 0.405(5)J. We find that even in the iTC phase, the superfluid order still persists (right inset of Fig. 4).

V. CONCLUSION AND OUTLOOK
We introduce a class of quantum many-body systems with oscillating retarded interactions in iTime and performed numerically exact QMC simulations. Our study finds an iTC phase with exotic zero-temperature properties and abnormal thermal behavior. To conclude this paper, we wish to outline some differences between our results and other related studies recently. First, our model is conveniently presented by a Euclidean action, so it is beyond the scope of Hamiltonian problems and does not contradict the "no time crystal ground state" theorem [17,18]. Second, even though it corresponds to a "time"-dependent problem, the "time" dependence enters our model through the relative time difference in terms of two-body interactions, in contrast to the early Floquet models whose time-dependence is explicitly through overall time [6,8]. This difference fundamentally changes the symmetry and a continuous symmetry model opens the wonder of what the analogue of Goldstone modes in spatial crystals, i.e, the ubiquitous and all important phonons, will be after the spontaneous breakdown of time translational symmetry. The potential of discovering "time" phonons is now open for further work. Finally, a thermodynamic physical system is found to support our "theoretical" time-crystal model. Standard statistical physics is applied to show the correspondence of our imaginary-time quantum mechanical model with retarded interaction (dependent on relative time!) to the ensemble of hard-core bosons embedded in engineered baths. This report therefore extends the study of time crystal to the Euclidean spacetime corresponding to thermal ensemble systems.
Our model can be derived from open quantum systems that have been studied extensively, such as quantum optomechanical resonators, QED cavities of interacting light and atoms, trapped ions [29] and Rydberg atoms [30]) coupled to synthetic thermal baths, etc. The crucial ingredient to realize the proposed idea is to look for a system that has retarded interaction and allows the system-bath coupling to be easily tuned. The above main conclusions are firmly supported by numerical exact QMC simulations of bosonic models with oscillating retarded interactions in relative iTime. Our method can be straightforwardly generalized to other bosonic and spin systems with retarded interactions, in one or higher dimensions. The competition between time-oscillating interaction and quantum-thermal fluctuations is expected to give rise to other forms of nontrivial temporal-thermal orderings that are completely absent in the conventional Hamiltonian systems. A generalization to fermionic models may be even more interesting since it opens up new possibilities of finding novel "non-fermi liquid". It is also interesting to study the real time counterpart of our model, which may be derived by integrating out a nonequilibrium bath. This may provide a new perspective for current time crystal researches since such an effective action has the full continuous (rather than discrete) symmetry in temporal translation. to a particle picture, i.e., the gas of instantons and antiinstantons (also referred in literature as to kinks and antikinks) with the same interacting energy [37]: where m(τ i ) = ±1 indicates a kink or antikink with the position τ i . The interacting energy V (τ − τ ′ ) between a pair of kinks can be directly evaluated by performing integral over τ and τ ′ , which gives rise to V (τ i − τ j ) = αe − (|τ i −τ j |−τ 0 ) 2 2σ 2 . (B3) Obviously, for those point-like objects (kinks and antikinks), this kind of interaction favor a crystalline pattern in τ -direction with lattice constant τ 0 , which corresponds to the iTC pattern in the world-line configurations.
To incorporate the quantum fluctuations in this picture, we perform QMC simulations of the model that is similar to what is discussed in the main text (hop-ping+retarded interaction), the only difference being that the kernel function replaced by Eq. (B1). In Fig. 5(c), we plot the onsite temporal density-density correlations, which shows a long-range crystalline order in τ -direction for a strong retarded interaction (α = 6J), indicating a spontaneous translational symmetry breaking in imaginary time. The lattice constant of the iTC is determined by the characteristic itime length (the minimum position) of the retarded interaction (τ 0 ). This example indicates that the oscillating feature of the retarded interaction is not necessary for the iTC phase, and hence proves that the phase can be observed in much boarder situations.