Divergence of thermalization rates driven by the competition
between finite temperature and quantum coherence

With the ramping of an external parameter, a classical isolated system evolves and attempts to thermalize into an equilibrium state at each specific number of the external parameter [1]. If the ramping is very slow, the system will closely follow the ramping and remain in a steady state. However, as the ramping gets faster, the internal thermalization cannot keep pace with the change of the external parameter, resulting in a time delay of the observable values reaching to the corresponded one at the steady state under each external parameter.

The scenario of a quantum isolated system at zero temperature is different. The phenomenon is fully described by the coherent evolution of wave functions under quantum mechanics [2, 3, 4], which is a different framework compared to thermodynamics. Nevertheless, as the particle number or system size increases, most systems can still exhibit thermalization described by the concept of entanglement thermalization [5, 6, 7, 8]. In this framework, statistical physics can provide consistent predictions as quantum mechanics at zero temperature [9, 10, 11, 12, 13, 14].

However, there remains a disparity between these two systems. In the case of a quantum system with finite temperature, the thermal effect interacts with and competes against quantum effects, but there have been few investigations in this context [15, 16, 17]. The quantum properties prevent the direct application of classical thermodynamics, and the finite temperature poses challenges in accurately calculating the quantum properties of these systems. This motivates us to use the ultracold-atom platform to investigate the thermalization behaviors under finite temperature, which provides a direct competition between quantum coherence and thermodynamics. In our study, we note thermalization rate divergence near absolute zero temperature caused by boson stimulation, in contrast to the convergence triggered by finite temperature and many-body interactions.

Our system is based on the Bose-Hubbard model with ultracold atoms [18, 19]. In Fig. 1a, we present the phase diagram of the homogeneous Bose-Hubbard model as a function of the temperature $T$ and the ratio $g=J/U$ between tunneling amplitudes $J$ and the on-site interaction strength $U$. The value of $g$ can be adjusted by manipulating the trap depth of optical lattices [20]. At zero temperature, the system exhibits a quantum phase transition separating the superfluid and Mott insulator. When the temperature is non-zero, the boundary of the superfluid shrinks which separates the superfluid and normal fluid with a continuous crossover [21]. The phase structure of superfluid/normal fluid bears similarities to the superconductor-metal transition.

Our measurement focuses on the ramping dynamics to investigate how the relaxation time $\tau_{SF}$ changes with the ramping speed $k$ at one particular parameter point $g_{0}$ under the influence of finite temperature effects. Inspired by the scaling relation in the Kibble-Zurek mechanism [22, 23, 24], we extend the concept of critical exponents to thermalization rate, characterizing the system’s thermalization property and how it responds to external ramping.

In our experiment, we non-adiabatically change the parameter $g$ with a ramping speed $k$, and measure the relaxation time at selected point $g_{0}$. When the system exhibits fast thermalization, the adiabatic theorem applies, and atoms adiabatically follow the changes in $g$, resulting in the same relaxation time $\tau_{SF}$ as the time required to ramp the parameter, that is, $\tau_{SF}=|g_{0}-g_{i}|/k$ where $g_{i}$ is the initial value of the ramping parameter. Consequently, $\mathcal{C}\rightarrow\infty$ approaches infinity for rapid thermalization. On the other hand, in the case of a system with a very slow thermalization rate, the observables do not track the external parameter ramping, and the relaxation time remains large but independent of ramping speed $k$, which leads to $\mathcal{C}\rightarrow 0$. This approach allows for a direct comparison of the intrinsic thermalization behaviors among different quantum phases with various $g_{0}$.

We utilize the improved band mapping technique [25, 26, 27, 28] to measure the quasi-momentum distribution of the system. This allows us to determine the number of atoms that exhibit coherence within the first Brillouin zone. This distribution directly corresponds to the phase coherence associated with superfluidity, which serves as an order parameter distinguishing superfluid and incoherent phases. We define the relaxation time $\tau_{SF}$ as the duration required to attain the same coherence level in an adiabatic case, and obtain the thermalization rate $\mathcal{C}$. This process is repeated for atoms with the same atom number $N_{0}$ but different initial temperatures $T_{0}$.

We prepare a total of $N_{0}=1.3(1)\times 10^{5}$ rubidium-87 atoms in condensates within a dipole trap, where the initial temperature $T_{0}$ can be adjusted. Here $T_{0}$ is measured by characterizing the thermal component of condensates (See Supplementary materials Section 2.). The dipole trap is a crossed dipole trap with the lowest vibrational frequencies of 42 and 60 Hz along horizontal and vertical directions. Then we adiabatically transfer the atoms from the dipole trap to the three-dimensional optical lattices, where the optical lattices are formed by the interference of three retro-reflecting laser beams with a wavelength at $\lambda=1064$ nm. The dipole trap is fully turned off when the optical lattices are turned on, and the trap depth of lattices is set at 5$E_{r}$ where $E_{r}=h^{2}/2m\lambda^{2}=h\times 2$ kHz is the recoil energy. This transfer process takes approximately 80 ms. Once loaded into the lattices, the atoms have a radius of around 10 $\mu$m. Subsequently, we hold the atoms in the lattices for 100 ms. The Gaussian beam waist of each lattice beam is 150 $\mu$m. As a result, the system is described by a Bose-Hubbard model with a harmonic-trap background, as [29, 30]

Here $a_{i}$ represents the annihilation operator of one particle at lattice site $i$, and $n_{i}=a^{\dagger}_{i}a_{i}$ is the particle number operator. The Gaussian shape of the lattice beams provides an isotropic harmonic trap, and the vibrational frequency $\omega$ is $2\pi\times$44 Hz at the trap depth $V=5E_{r}$. The frequency $\omega$ scales as $\omega\propto\sqrt{V}$. Due to the Gaussian envelope of the lattice beams, the variation of $J/U$ for atoms on the outer regions compared to the center part is only 1% difference [31, 32].

Firstly, we calibrate our system (Fig. 2a) [25, 26, 27, 28] to distinguish the different phase components and determine the number of incoherent atoms. We employ the relative incoherent fraction $\gamma_{inc}$ to characterize the system. Since the direct measurement of the real temperature in optical lattices is challenging [33], we use the initial temperature $T_{0}$ as a characterization parameter instead. In Fig. 2b, we present the data of band mapping as a function of $T_{0}$ at the trap depth $V=5E_{r}$. As $T_{0}$ increases, the thermal effects become more dominant, resulting in a higher fraction of incoherent atoms. However, when $T_{0}$ exceeds 220 nK, atoms are excited into the second Brillouin zone, where the incoherent fraction becomes difficult to extract, which is what we want to avoid. At around 50 nK, almost all the atoms remain coherent and stay in the superfluid phase at a trap depth of $5E_{r}$. In addition, we also calibrate the change of the incoherent fraction $\gamma_{inc}(T_{0})$ at different points $g_{0}$ corresponding to different trap depth $V$ (Fig. 2c and d). Here when $V$ is above $13E_{r}$, the Mott insulators start to gradually appear for Rubidium-87 atoms [20, 32] in such inhomogeneous optical lattices. Based on these calibrations, we obtain the steady incoherent fraction $\gamma_{inc}(T_{0})$ versus the initial temperature $T_{0}$ and different trap depth $V$.

Next, we perform the ramping of the trap depth $V$ linearly with a ramping speed $k$ for samples with different initial temperatures $T_{0}$ (Fig. 3a). At each specific $V$, we perform the band mapping and extract $\gamma_{inc}(V,T_{0},k)$ (Fig. 3b to e). During this process, we ensure that there are no excitations into the second Brillouin zone when controlling $T_{0}$ below 220 nK. When the temperature is low (45 nK in Fig. 3b), there is a clear turning point at 13$E_{r}$ for different ramping speed $k$. The dynamical responses are almost the same before the critical point $13E_{r}$ and start to furcate after it. Therefore, we utilize the steady incoherent fraction $\gamma_{inc}(T_{0})$ at the critical point $13E_{r}$ as a reference to measure the time required to reach the same level of incoherence. This procedure provides us with a relation time $\tau_{SF}(T_{0},k)$ for each ramping speed $k$ and $T_{0}$. By plotting the logarithmic relation of $\tau_{SF}(T_{0},k)$ versus $k$ in Fig. 3f and fitting the data using the scaling relation $\tau_{SF}(T_{0},k)\propto k^{-\alpha}$, we determine the value of $\alpha$, which is related to the thermalization rate $\mathcal{C}$ through $\mathcal{C}=\alpha/(1-\alpha)$. particularly, the critical exponent $\nu z$ can be obtained via $\nu z=\mathcal{C}$ at the critical point. In the logarithmic plot, the fitting lines exhibit good agreement with the data, confirming the validity of the scaling relation. Specifically, at $T_{0}=45$ nK, $\alpha$ is 0.92(1) corresponding to $\mathcal{C}=12.05^{{}_{+1.99}}_{{}^{-1.52}}$. At $T_{0}=76$ nK, $\alpha$ is 0.85(2) and the corresponded $\mathcal{C}$ is $5.72^{{}_{+0.92}}_{{}^{-0.72}}$. At $T_{0}=128$ nK, $\alpha$ is 0.71(3) and the corresponded $\mathcal{C}$ is $2.42^{{}_{+0.45}}_{{}^{-0.35}}$. In Fig. 3g, we plot the relation between $\mathcal{C}$ and $T_{0}$, while the inset panel shows the plot of $\alpha$ versus $T_{0}$. It is evident that $\mathcal{C}$ increases dramatically as the temperature approaches zero. We fit this change using a trial function, which only serves as an empirical fit based on the data

We think the superfluid state itself should bear the responsibility of divergence instead of any parameters in Hamiltonian since only incoherence can converge the thermalization rates. In superfluid, most of the atoms condense into the zero-momentum state, and this leads to the bosonic stimulation [34, 35] that the excitation speed from superfluid will be amplified by the condensed atoms with a factor $\sqrt{\tilde{N}}$. Here $\tilde{N}$ is the condensed-atom number. When the temperature increases, the increasing proportion of incoherence attenuates the bosonic stimulations and quantum coherence, and then the thermalization rates converge.

One way to further support this argument is to inspect the thermalization rates of Mott insulators. Even at zero temperature, the superfluid will diminish when the interaction dominates the system, and this will also suppress the bosonic stimulation. Therefore, we examine the thermalization rate $\mathcal{C}$ using different truncation trap depth. In Fig. 4, we illustrate the variation of $\mathcal{C}$ as a function of trap depth for different temperatures. We observe that the divergence of $\mathcal{C}$ exhibits a similar pattern to the phase diagram. Specifically, the superfluid phase is associated with larger $\mathcal{C}$ in the regime with lower $T_{0}$ and smaller $V$. In the Mott insulator or normal fluid regime, $\mathcal{C}$ rapidly decays to around 2. This suggests that the divergence of $\mathcal{C}$ is an intrinsic behavior of superfluid due to the quantum coherence instead of disorders.

In conclusion, our study examines the thermalization behaviors of cold bosonic systems in optical lattices at various temperatures. Based on the scaling relations of relaxation time versus the external parameter ramping rate, we extract the thermalization rate which exhibits a divergence as the temperature approaches zero, reflecting the intrinsic stimulation of the superfluid phase. Conversely, the thermalization rate converges as either the temperature or interaction increases, where the phase decoherence or thermal excitation suppresses the speeding up. This is a direct competition between quantum and thermal effects. We believe our work provides fresh perspectives on quantum systems at finite temperature and related to the divergent [36, 37, 38, 39] or scaling [40, 41, 42, 43, 44, 45] behaviors in other quantum materials.

We acknowledge the stimulating discussions from Zhiyuan Yao. This work is supported by National Key Research and Development Program of China (2021YFA0718303, 2021YFA1400904), National Natural Science Foundation of China (92165203, 61975092, 11974202), and Tsinghua University Initiative Scientific Research Program.

In the experiment, we use the improved band-mapping method to measure the incoherent fraction $\gamma_{inc}$. The band-mapping method is an efficient approach to measuring the quasi-momentum distribution for different phases. By adiabatically ramping down the power of all lattice beams, the lattice quasi-momentum is converted to real-space momentum. Nonetheless, during this process, the momentum distribution is disturbed by the on-site interaction, thus we cannot obtain the ideal uniform distribution for Mott insulators. In our improved method, unlike the slowly ramped $x$ and $y$ lattices, the $z$ lattice is turned off immediately to release the interaction, and this gives us a better converted real-space momentum.

In the deep Mott region, the wavefunction of each atom is a Wannier function localized within one site and is also a superposition of all possible Bloch functions within one Brillouin zone. Therefore, the quasi-momentum is uniformly distributed as a flat plateau in the first Brillouin zone. On the contrary, the superfluid phase corresponds to a central peak since the atoms tend to be in the zero-momentum state. We define the incoherent fraction as:

where $A_{pk}$ denotes the area of the central coherent peak, and $A_{tot}$ is the total area of the first Brillouin zone. In our experiment, the phase decoherence leads to the increase of incoherent atoms, thus the increase of $\gamma_{inc}$.

In order to eliminate the impact of in-alignment errors, the quasi-momentum distribution is firstly centralized and symmetrized. We choose an approximate central momentum and turn the right momentum distribution to the left, then calculate the residual error in a chosen area with the unprocessed data. The central point is swept in a reasonable range to find the optimal point where the residuals are minimal. After the centralization, the coherent peak is easy to be distinguished based on the plateau.

In this section, we show how we measure the initial condensate temperature $T_{0}$. To measure $T_{0}$, we release the condensate from the optical dipole trap for a 30-ms-long free expansion, and then apply a resonant absorption imaging. Based on the time-of-flight method, we obtain the two-dimensional density $n(x,y)$. Then, we integrate one direction and make the one-dimensional density $n(x)$ to be fitted by a bimodal distribution [46, 47],

where $x_{0}$ stands for the central place of atoms, $n_{th}^{0}$ and $n_{BEC}^{0}$ are the central densities of thermal atoms and condensates respectively, $\sigma_{th}$, $\sigma_{BEC}$ are the radii of thermal atoms and condensates, and $g_{\alpha}(x)=\sum_{n}z^{n}/n^{\alpha}$ is the Bose function.

Then we calculate the condensate fraction $N_{0}/N$ using the fitting results above and the critical phase transition temperature $T_{c}$ by the formula

where $\omega$ is the vibrational frequency of the trap, $N$ is the total atom number, and $k_{B}$ is the Boltzmann constant. According to the relation of the condensate fraction and the temperature

where $a$ equals to 3 for the harmonic trap, we obtain the temperature corresponding to each experimental measurement.

The fitting process consists of two steps. In the first step, the data of $\gamma_{inc}$ under the same temperature is interpolated into smooth curves. Then we take $\gamma_{inc}$ at the critical point as a criterion to extract the relaxation time. The details are described in the following subsections.