Imaginary Time Crystal of Thermal Quantum Matter^*

Introduction.—The phases of quantum matter are, in general, categorized according to symmetries and their spontaneous breaking. A striking possibility was put forward,^[1] dubbed "time crystal", where the interacting particles spontaneously organize themselves into a "time" coherent state that breaks temporal translational symmetry. Temporal periodicities spontaneously emerge analogous to spatially ordered crystals. This phase has raised considerable interest,^[2–17] and was subsequently proven to be forbidden in thermodynamic equilibrium states.^[18,19] A variant of the phase was proposed^[6–8] by considering certain non-equilibrium condition, and has quickly been observed in experiments.^[20,21]

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, finite-temperature 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 exists and what novel effects should manifest in ground state and thermal response functions.

In this Letter, 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. The feature satisfies the two conditions that may previously seem to be 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.

Model and method.—Our model is hard-core bosons in either one-dimensional (1D) lattice of L sites with system Hamiltonian: ${H}_{{\rm{s}}}=\displaystyle \sum _{\langle i\rangle }-J[{a}_{i}^{\dagger }{a}_{i+1}+h.c.]$, where J denotes the hopping amplitude. Focus is on the half-filling case throughout this study. On each site i the density operator of the hard-core boson n_i additionally couples to a local bath. We assume that 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^−βH_tot (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 (see Supplemental Material (SM) for details)

$$\begin{eqnarray}&&{\hat{\rho }}_{{\rm{s}}}=\displaystyle \frac{1}{{Z}_{{\rm{s}}}}{e}^{-\beta {H}_{{\rm{s}}}-{S}_{{\rm{R}}}},\end{eqnarray} \tag{ 1 }$$

where H_s is the system Hamiltonian and ${Z}_{{\rm{s}}}={{\rm{Tr}}}_{{\rm{s}}}{\hat{\rho }}_{{\rm{s}}}$. S_R describes the effective action form of the bath induced retarderd interaction which is onsite in space and translationally invariant in iTime:

$$\begin{eqnarray}&&{S}_{{\rm{R}}}=-\displaystyle {\int }_{0}^{\beta }{\rm{}}d\tau \displaystyle {\int }_{0}^{\beta }{\rm{}}d{\tau }^{^{\prime} }\displaystyle \sum _{i}[{n}_{i}(\tau )-{n}_{0}]D(\tau -{\tau }^{^{\prime} })[{n}_{i}({\tau }^{^{\prime} })-{n}_{0}],\end{eqnarray} \tag{ 2 }$$

where n₀ = 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 ω₀ gives rise to the kernel function D(τ − τ') ∼ e^{−ω₀|τ − τ'|} + e^{−ω₀(β −|τ − τ'|)} for large β. We will show later that the iTC phase can exist for a broad class of kernel functions. Here we choose

$$\begin{eqnarray}\begin{array}{lll}D(\tau -{\tau }^{^{\prime} }) & = & \alpha [F(\tau -{\tau }^{^{\prime} })+F(\beta -|\tau -{\tau }^{^{\prime} }|)],\\ F(\tau -{\tau }^{^{\prime} }) & = & {e}^{-{\omega }_{d}|\tau -{\tau }^{^{\prime} }|}\cos [2\pi {\omega }_{c}(\tau -{\tau }^{^{\prime} })],\end{array}\end{eqnarray} \tag{ 3 }$$

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, a bosonic bath induced interaction 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^[22–24]). Now we need an RKKY-like interaction, but with oscillating decay in iTime instead of real space. Consider

$$\begin{eqnarray*}&&D(\tau )=\displaystyle \sum _{m}\displaystyle \frac{{e}^{i{\omega }_{m}\tau }}{i{\omega }_{m}-{\Sigma }_{b}(i{\omega }_{m})}\end{eqnarray*}$$

with ω_m = 2π m/β 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,^[25] in the first section of SM, we propose a similar microscopic bath model that can induce the oscillating retarded interaction as expressed 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 the second section of SM, 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. [26], we generalize the QMC simulation with worm-type updates^[27] 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 interaction [Eq. (2)] only appears in the diagonal parts (under the Fock basis) of the effective action, they do not cause extra minus sign problem. Thus our QMC simulations can be considered to be numerically exact since our "system" Hamiltonian (of hard-core bosons) is also positive definite.

Groundstate properties.—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, 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 quasi-superfluid 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 interaction leads to novel quantum phases.

The crystalline order can be identified from the density correlation functions in space and iTime:

$$\begin{eqnarray*}&&S(r,\tau )=\displaystyle \frac{1}{L}\displaystyle \sum _{i}\langle ({n}_{i}(0)-{n}_{0})({n}_{i+r}(\tau )-{n}_{0})\rangle .\end{eqnarray*}$$

We first focus on the onsite correlations (r = 0). As shown in Fig. 1(a), for small α, S(0,τ) still decays algebraically with τ, but is accompanied by oscillations. For 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}_{\tau }^{2}=\displaystyle \frac{1}{\beta }\displaystyle {\int }_{0}^{\beta }d\tau S(0,\tau )\cos [2\pi {\omega }_{c}\tau ]$. A finite-size scaling of m_τ is shown in Fig. 1(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 phase transition can be 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,τ)). The normalized correlation length ξ_τ/β as a function of α for different system sizes and β is plotted in Fig. 1(e), where we find a crossing point at ${\alpha }_{c}^{\tau }=0.23(1)$, indicating a scaling invariant critical point. In the iTC, ξ_τ/β linearly scales with the system size, and extrapolates to a finite value in the thermodynamic limit (inset of Fig. 1(e)).

Zoom In Zoom Out Reset image size

For sufficiently large α, it is interesting to notice that the crystalline pattern emerges not only in iTime, but also in real space, which can be identified by the equal-time correlation function S(r,0) as well as its structure factor $S(Q)=\displaystyle \frac{1}{L}\displaystyle \sum _{r}{e}^{iQr}S(r,0)$. In the case of half-filling, the spatial crystalline pattern is a charge-density-wave (CDW) state with a periodicity twice that of the lattice, as shown in Fig. 1(b). A finite size scaling of the order parameter ${m}_{{\rm{s}}}=\sqrt{S(Q=\pi )}$ is shown in Fig. 1(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 hard-core bosons being separated as far as possible to avoid blocking the oscillations of each other. Therefore at half-filling the tendency towards CDW order is expected in real space. Similar to the above analysis, we also plot the normalized correlation length ξ_s/L in real space as a function of α in Fig. 1(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,\tau )=\langle {a}_{i}^{\dagger }(\tau ){a}_{i+r}(0)\rangle$. The equal-time part G(r,0) is shown in Fig. 2(a), which decays algebraically with r (accompanied by an oscillation) even for large α, indicating a quasi-superfluid order in the iTC. This order can be further verified from the superfluid density, which can be calculated as ${\rho }_{{\rm{s}}}=\displaystyle \frac{L}{\beta }\langle {W}^{2}\rangle$, 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. 2(a), ρ_s grows monotonically with α, indicating that the retarded interaction enhances the mobility of bosons. Now we turn to the temporal part of G(r,τ), and focus on its k = 0 component: $G(k=0,\tau )=\displaystyle \frac{1}{L}\displaystyle \sum _{r}G(r,\tau )$. Figure 2(b) shows that G(k = 0,τ) decays exponentially with τ, indicating a finite single-particle charge gap. This agrees with the results of the compressibility $\kappa =\displaystyle \frac{\beta }{L}(\langle {N}^{2}\rangle -{\langle N\rangle }^{2})$ with N the total particle number. Figure 2(c) shows that κ decays to zero for large α, thus the iTC is incompressible.

Zoom In Zoom Out Reset image size

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 possesses 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.^[28,29] Since a world line needs to be continuous, the temporal-spatial configurations minimizing the retarded interacting action (2) are highly degenerate. As shown in Fig. 2(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. 2(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. 2(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.

Temperature crystal.—The most striking effect of iTC 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 behavior 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. 2(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.

Experimental realization.—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. Let us elaborate on what was expressed when introducing the model Hamiltonian. Such an interaction can be induced by a fermionic bath with finite life time, which generally arises from fermion self-scattering or from the scattering of bath fermions with impurities or phonons. Instead of the non-Hermitian mechanism as proposed in Ref. [25] (where the effective interaction is also induced by the scattering with phonons albeit they had a completely different purpose), here we propose a strongly interacting ultracold fermionic bath, whose self-energy naturally contains an imaginary part (finite line width in spectrum). Recent breakthrough in an ultracold atomic system^[30] has shed light on realization of such an unconventional interaction. In Ref. [30], fermion-mediated interactions have been observed in a bose-fermi mixture in the regime that the dynamics of fermions is much faster than that of the bosons. The difference needed for our model appears to adjust the characteristic time scales of the bosons and fermions to be compatible, so that the time-delayed feature of the bath induced interaction can not be neglected.

Conclusion and outlook.—We have introduced 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 study, we wish to outline some differences between our results and other related studies recently. Firstly, 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.^[18,19] Secondly, 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 the 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 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 non-equilibrium 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.

Acknowledgment