Quantum Thermalization and the Expansion of Atomic Clouds

TL;DR

This paper investigates whether quantum many-body systems can exhibit behavior fundamentally different from classical systems, focusing on atomic cloud expansion and identifying conditions where quantum effects induce diffusive dynamics.

Contribution

It demonstrates that quantum correlations near Bose-Einstein condensation can alter atomic cloud expansion, showing a transition from ballistic to diffusive behavior.

Findings

Quantum correlations influence cloud expansion near condensation temperature.

Most quantum dynamics are indistinguishable from classical ballistic motion.

A notable quantum effect causes diffusive behavior in atomic clouds.

Abstract

The ultimate consequence of quantum many-body physics is that even the air we breathe is governed by strictly unitary time evolution. The reason that we perceive it nonetheless as a completely classical high temperature gas is due to the incapacity of our measurement machines to keep track of the dense many-body entanglement of the gas molecules. The question thus arises whether there are instances where the quantum time evolution of a macroscopic system is qualitatively different from the equivalent classical system? Here we study this question through the expansion of noninteracting atomic clouds. While in many cases the full quantum dynamics is indeed indistinguishable from classical ballistic motion, we do find a notable exception. The subtle quantum correlations in a Bose gas approaching the condensation temperature appear to affect the expansion of the cloud, as if the system has turned into a diffusive collision-full classical system.

Equations

$$\mathcal{M}=-\log \rho .$$

$${\rho }_0=\frac{1}{Z_A{Z}_B}{e}^{-{\beta }_A{\mathcal{H}}_A}\otimes e^{-{\beta }_B{\mathcal{H}}_B}$$

$$\mathcal{M}(t)=e^{-i\mathcal{H}t}{\mathcal{M}}_0{e}^{i\mathcal{H}t}.$$

$${\mathcal{M}}_0=\sum _{k{k}^{\prime }}{m}_{k{k}^{\prime }}{\hat{c}}_k^{†}{\hat{c}}_{k^{\prime }}^{}+\log Z.$$

$${\hat{c}}_k^{†}(t)=e^{-i\mathcal{H}t}{\hat{c}}_k^{†}{e}^{i\mathcal{H}t}=e^{-i{\xi }_{k}t}{\hat{c}}_k^{†},$$

$$\mathcal{M}(t)$$

$$\equiv$$

$$\hat{G}(t)={\left[e^{\hat{m}(t)}-\eta \right]}^{-1}.$$

$$\mathcal{H}=-J\sum _{\mathbf{r},\delta }\left({\hat{c}}_{\mathbf{r}+\delta }^{†}{\hat{c}}_{\mathbf{r}}^{}+{\hat{c}}_{\mathbf{r}}^{†}{\hat{c}}_{\mathbf{r}+\delta }^{}\right)=\sum _{\mathbf{k}}{ϵ}_{\mathbf{k}}{\hat{n}}_{\mathbf{k}}$$

$${\mathcal{H}}_A=-J\sum _{r_x,r_y,\dots =1}^{L_A-1}\sum _{\delta }\left({\hat{c}}_{\mathbf{r}+\delta }^{†}{\hat{c}}_{\mathbf{r}}^{}+{\hat{c}}_{\mathbf{r}}^{†}{\hat{c}}_{\mathbf{r}+\delta }^{}\right).$$

$$m_{j\ell }(t)={\beta }_B{\delta }_{|j-\ell |=1}+({\beta }_A-{\beta }_B)\int \frac{d^d\mathbf{k}{d}^d{\mathbf{k}}^{\prime }}{(2\pi {)}^{2d}}{e}^{-i\mathbf{k}{\mathbf{r}}_j+i{\mathbf{k}}^{\prime }{\mathbf{r}}_{\ell }}\left(\sum _{i=1}^d{e}^{i{k}_i}+e^{-i{k}_i^{\prime }}\right)\left(\prod _{i=1}^d\frac{e^{i(k_i-k_i^{\prime })(L_A-1)}-1}{e^{i(k_i-k_i^{\prime })}-1}\right)e^{i({ϵ}_{\mathbf{k}}-{ϵ}_{{\mathbf{k}}^{\prime }})t}.$$

$$\Delta m_{j,j+1}(t\gg 1)=2\Delta \beta (0){\left(\frac{L_A-1}{2\pi Jt}\right)}^d\sim \frac{V_A}{t^d}.$$

$$j_E(x)=-\kappa (T(x))\nabla T(x)$$

$${\partial }_{t}T=\mathcal{D}{\nabla }^{2}T.$$

$$\Delta T(x,t)=\frac{1}{2}(T_A-T_B)\left(\mathrm{Erf}\left[\frac{a-2x}{4\sqrt{\mathcal{D}t}}\right]+\mathrm{Erf}\left[\frac{a+2x}{4\sqrt{\mathcal{D}t}}\right]\right).$$

$$\Delta T(x=0,t\gg 1)\sim (T_A-T_B){\left(\frac{a}{2\sqrt{\pi \mathcal{D}t}}\right)}^d\sim \frac{V_A}{t^{d/2}}.$$

$$n(x,v,t)=n(x-vt,v,0).$$

$$n(x,v,0)=\frac{1}{2\pi }\frac{dk}{dv}\frac{1}{e^{\beta (x)(ϵ(v)-\mu (x))}-1}$$

$$n(x=0,t)$$

$$\approx$$

$$S=-\mathrm{Tr}\rho (t)\log \rho (t)=\mathrm{Tr}\mathcal{M}(t)\rho (t).$$

$$S=\sum _{\alpha }\left[-g_{\alpha }\log g_{\alpha }+\eta (1+\eta g_{\alpha })\log (1+\eta g_{\alpha })\right],$$

$$\mathcal{H}=\int dx\left({\psi }_R^{†}(x)iv{\partial }_x{\psi }_R(x)-{\psi }_L^{†}(x)iv{\partial }_x{\psi }_L(x)\right)$$

$${\psi }_R^{†}(x,t)=\int \frac{dk}{2\pi }{\psi }_{R,k}^{†}{e}^{ik(x-vt)}.$$

$${\mathcal{M}}_0={\beta }_B\mathcal{H}+({\beta }_A-{\beta }_B){\int }_0^{L_A}dxh(x)$$

$${\int }_0^{L_A}dx{h}_R(x,t)={\int }_{vt}^{L_A+vt}dx{\psi }_R^{†}(x)iv{\partial }_x{\psi }_R(x)$$

Full text

Introduction

The laws describing classical gases, most notably the Second Law of Thermodynamics, seem at odds with the principle of unitary time evolution in quantum physics.[1] However, high energy states are densely many-body entangled and consequently the Eigenstate thermalization hypothesis (ETH)[2, 3, 4] claims that the outcomes of local measurements will be at long times indistinguishable from the outcome of the measurement in a thermal mixed state, at a temperature consistent with the energy of that state.[5, 6]

Is this also true for a cloud of non-interacting quantum particles confined in a potential, which is suddenly released and allowed to expand in an infinite bath? This is actually similar to the key ‘time-of-flight measurement’ in many cold atom experiments.[7, 8] After suddenly releasing the confining potential the atomic clouds expand, and by assuming that this is governed by ballistic, collision-less atomic motion the initial velocity distributions can be deduced from the expansion of the cloud. Invariably, it has been assumed that this expansion is governed by a purely classical Newtonian or wave kinematics, and this is undoubtedly a correct procedure to follow.

However, it is not at all obvious why this works. After all, before releasing the trapping potential, one may be in a quantum regime with Bose condensation or Fermi-degeneracy. How can these atoms suddenly behave like classical canon balls? In the next section, we will present a method to compute local observables exactly in the full quantum evolution by evaluating the logarithm of the density matrix. A first result is shown in Fig.1: under the conditions of the cold atom experiments the full quantum dynamics is indeed indistinguishable from classical ballistic expansion.

We then address cooling, where the atoms are released in a particle bath which is at a lower temperature than the trapped particles. When the temperature of the bath is high enough we find an expansion consistent with the classical expectation: since the particles do not collide, the hot cloud cools ballistically. Similarly, when the cloud and the bath are both formed from fermions the system behaves classical. However, for a cloud of bosons cooling into a bosonic bath at a temperature approaching the condensation temperature, the cooling is governed by diffusion! In Fig. 2 we show how the energy density of a ‘hot’ cloud in a cold bath spreads out in time, marking a clear difference between classical diffusion, ballistic fermionic behavior and again diffusion for a low temperature bosonic bath. Quantitatively, the difference between ballistic and diffusive behavior can be shown by measuring the total energy density $\Delta E$ in the region of the original cloud relative to the bath energy density, as shown in Fig. 3: ballistic decay is characterized by $\Delta E\sim t^{-d}$ whereas diffusion satisfies $\Delta E\sim t^{-d/2}$.

This is our main result. We have identified a circumstance where the quantum evolution becomes sharply distinguishable from the analogous classical evolution. In the classical system diffusional expansion requires collisions, but these are collision-less quantum particles. We will explain how to test this prediction in cold atom experiments, but first we elucidate how these matters are computed.

Method

The traditional approach to evaluate quantum time evolution is by repeated application of the time evolution operator $e^{-i\mathcal{H}dt}$ with small temporal steps $dt$. However, with this procedure it is impossible to reach times later than $t\sim 1/E$, where $E$ is a typical energy scale of the system. The hypothesis of thermalization provides us now with a simpler way to compute time evolution through the modular Hamiltonian $\mathcal{M}$, which is the logarithm of the density matrix,

[TABLE]

As we will see, at late times $\mathcal{M}$ will simplify dramatically. Since we are interested in a hot cloud in a cold bath, our initial density matrix will have the form111This is equivalent, up to boundary terms, to $\rho_{0}\sim\mathrm{Tr}_{B}\,e^{-\beta_{A}\mathcal{H}}\otimes\mathrm{Tr}_{A}\,e^{-\beta_{B}\mathcal{H}}$.

[TABLE]

where $\mathcal{H}_{X}$ and $\beta_{X}$ are the total Hamiltonian and inverse temperatures respectively, restricted to the subsystems $X=A,B$. The time evolution of the modular Hamiltonian follows directly from the von Neumann equation for the time evolution of the density matrix,

[TABLE]

For noninteracting systems $\mathcal{H}=\sum_{k}\xi_{k}n_{k}$, the initial modular Hamiltonian following from Eqn. (2) can be written as,

[TABLE]

The modular matrix $\hat{m}=m_{kk^{\prime}}$ is Hermitian; the sum runs over the momenta $k$ of the particles, while the constant $\log Z=-\eta\mathrm{Tr}\;\log\left(1-\eta e^{-\hat{m}}\right)$, with $\eta=-1$ for fermions and $\eta=+1$ for bosons. The time evolution of both fermion and boson field operators appearing in the modular Hamiltonian is for the free system simply given by,

[TABLE]

This implies for the time dependence of the modular Hamiltonian,

[TABLE]

It follows that time evolution corresponds with a unitary transformation on the modular matrix. Local observables as the occupation numbers and the energy are in turn functions of the equal-time Greens function at time $t$, $G_{ij}(t)=\mathrm{Tr}\;\hat{c}^{\dagger}_{i}\hat{c}^{\phantom{\dagger}}_{j}\rho(t)$, in terms of the modular matrix

[TABLE]

The advantage of this formulation starts to shimmer through. The intricacies of the full quantum evolution are absorbed in the strongly oscillating factors occurring in Eq. (7). These will rapidly average away such that in the limit $t\rightarrow\infty$, the modular Hamiltonian approaches the actual Hamiltonian, $\mathcal{M}(t)\rightarrow\beta\mathcal{H}$ when expressed in a local basis.

Expansion of a noninteracting hot gas in a cold bath

To see how this works let us consider some examples. Relativistic systems are discussed in the supplementary material[9], reproducing the wisdom that these thermalize instantaneously once full causal contact is established[1, 2, 3, 4]. To model non-relativistic atoms we resort to a lattice regularization in the form of a hypercubic lattice in $d$ dimensions with nearest neighbor hopping,

[TABLE]

with $\epsilon_{\bf k}=-2J\sum_{i=1}^{d}\cos k_{i}$. Given our initial hot cloud state the modular Hamiltonian equals $\mathcal{M}(t)=\beta_{B}\mathcal{H}+(\beta_{A}-\beta_{B})\mathcal{H}_{A}(t)$ where the Hamiltonian of the subsystem $A$ is at $t=0$ equal to

[TABLE]

Under time evolution this hot cloud spreads out and at $t>0$ we express $\mathcal{M}(t)=-J\sum_{j\ell}m_{j\ell}(t)\hat{c}^{\dagger}_{j}\hat{c}^{\phantom{\dagger}}_{\ell}$ in terms of the elements of the modular matrix $m_{j\ell}(t)$ in the real space basis,

[TABLE]

Recall that thermalization in the ETH sense implies that the second term should vanish at late times. Indeed, using the continuum approximation $\epsilon_{\bf k}\approx Jk^{2}-\mu$ for $t\gg 1/\epsilon_{\bf k}$, and thus $(\epsilon_{\bf k}-\epsilon_{{\bf k}^{\prime}})t\approx({\bf k}+{\bf k}^{\prime})({\bf k}-{\bf k}^{\prime})Jt$, we find for a site $j\in A$,

[TABLE]

Regardless of the statistics of the particles, the modular Hamiltonian approaches the final thermal state with a ballistic powerlaw decay.[9]

However, the experimentally relevant local energy density in subsystem $A$ can approach the bath value in different manners, pending the quantum statistics of the particles as illustrated in Fig. 3. Fermions are consistently subjected to a ballistic decay of the energy difference between the bath and the subsystem $A$ (Fig. 3a) and the resulting energy flow profile (Fig. 2c) displays a smoothened light-cone following the Lieb-Robinson bound with $v_{LR}=2J$.[14] Turning to bosons, the surprise we announced becomes manifest: we find a crossover from ballistic behavior at high bath temperatures to diffusive $\Delta E\sim t^{-d/2}$ at low bath temperatures. For both $d=1$ and $d=2$ dimensions (Fig. 3b,c), the crossover occurs around the point where the lattice thermal de Broglie wavelength corresponds to the interparticle spacing. This suggests that diffusive behavior is a consequence of the wave-like nature of the bosons, where $\partial_{t}\psi\sim\partial_{x}^{2}\psi$.

The energy profile of the diffusive case (Fig. 2d) is surprisingly reminiscent of the classical Fourier’s law of heat diffusion (Fig. 2a). However, one should not be fooled by this apparent relation to classical diffusion. After all, we are considering noninteracting particles and the equivalent classical description of our set-up is through a distribution of particles and velocities $n(x,v,t)$ that evolves ballistically $n(x,v,t)=n(x-vt,v,0)$. For the expansion into a cold bosonic bath the classical picture still yields a ballistic spread, while the exact quantum evolution displays diffusive behavior.[9] The diffusive behavior for cold bosonic baths is therefore a genuine quantum effect.

Conclusion and Outlook

This interesting crossover can be probed directly in experiments using cold atoms, following the protocol illustrated in Fig. 4.[15, 16, 17] Initially, one prepares a cloud of atoms tuned to be noninteracting using the Feshbach resonance. Using optical lattice techniques a barrier is created in between $A$ and $B$, and a separate laser excites $A$ to be at a different temperature than $B$. At time $t=0$ the barrier is removed and the system will evolve as described. To measure the energy density in subsystem $A$ after a time $t$, one reintroduces the barrier, let the atoms in the bath $B$ escape, followed by time-of-flight measurements of the distribution of momenta of the atoms in $A$. From the distribution of these momenta the total kinetic energy can be reconstructed. The experiment is then repeated to obtain the energy density in $A$ at every time instance. In this way the curves of Fig. 3, for either ballistic or diffusive behavior, can be experimentally measured.

It might be a surprise to observe thermalization in integrable non-interacting systems, but it is quite straightforward that this happens for local quenches such as the one studied here.[18, 19] Even though there are many integrals of motion, there is no conservation law that restricts certain degrees of freedom to remain within $A$. Note, however, that systems where the integrals of motion are truly local, as is the case for Anderson insulators[20] or the many-body localized phase[21, 22], information remains within $A$ and no thermalization will occur.

A critical reader might object that the system we study actually displays an entropy decrease. However, much like refrigerators, we reduce the entropy of subsystem $A$ by increasing the bath entropy by at least the same amount. In fact, while the total entropy remains constant in any quantum system, the mutual information $\mathcal{I}_{AB}(t)=S_{A}(t)+S_{B}(t)-S_{A\cup B}$ increases upon thermalization since the subsystem $A$ and the bath $B$ become entangled. This increase in mutual information should be considered the quantum version of the Second Law.[1] However, it remains an open question to prove this increase for thermodynamically large systems as the Second Law requires.

Acknowledgements

We are thankful to Tarun Grover, Laimei Nie, Mike Zaletel and Immanuel Bloch for discussions. L.R. was supported by the Dutch Science Foundation (NWO) through a Rubicon grant and by the National Science Foundation under Grant No. PHY11-25915 and Grant No. NSF-KITP-17-019.

Author contributions statement

L.R. performed the numerics, and L.R. and J.Z. wrote the manuscript together.

Additional information

The authors declare no competing financial interests.

Appendix A Thermalization of classical systems - Fourier’s Law

In the main manuscript we consider a hot system $A$ at temperature $T_{A}$ immersed in a cold bath at temperature $T_{B}$. How will thermal equilibrium be reached according to the classical theory of thermal diffusion? There it is assumed that a system is locally in thermal equilibrium, such that one can define a temperature $T(x)$ at each point in space. If there exists a temperature gradient, energy will flow from hot to cold according to Fourier’s law,

[TABLE]

where $\kappa$ is the thermal conductivity of the material and $j_{E}(x)$ is the energy current. If we are in a regime where both the specific heat $c_{V}$ and the thermal conductivity $\kappa$ are approximately independent of temperature, Fourier’s Law becomes a diffusion equation

[TABLE]

where the diffusion constant is $\mathcal{D}=\frac{\kappa}{c_{V}}$. This equation can be solved using the heat kernel.

Let us look explicitly at an initial state with a hot cloud at temperature $T_{A}$ for $|x|<a/2$, and and a bath at $T_{B}$ for $|x|>a/2$. The resulting solution of the heat diffusion equation yields

[TABLE]

In the left panel of Fig. 1 in the main manuscript, we show how the heat of our cloud spreads according to the above formula.

Note that the temperature difference at long times falls of in a power law fashion, $\Delta T(x=0,t\gg 1)\sim\frac{a(T_{A}-T_{B})}{2\sqrt{\pi\mathcal{D}t}}$. In higher dimensions $d$, the above equations straightforwardly generalize to

[TABLE]

Therefore, if the energy of temperature of a system decays as a powerlaw $t^{-d/2}$, we call this diffusion.

Appendix B Comparison between classical and quantum description

Now consider another classical system: a gas of collision-less non-interacting particles. At time $t=0$, we can characterize this gas as having a distribution of particles in position and velocity, $n(x,v,0)$. The particle density as a function of position is $n(x)=\int dv\;n(x,v,0)$, and with an energy per particle that depends only on velocity, $\epsilon(v)$, the energy density is given by $E(x)=\int dv\;\epsilon(v)n(x,v,0)$.

Because the particles are collision-less and have no further interactions, the velocity is conserved. This means that the full distribution at a later time $t$ can be expressed in terms of the initial distribution as

[TABLE]

To model a generic system of bosons as is done in the main manuscript, we can start with an initial distribution

[TABLE]

with $\epsilon_{\bf k}$ the boson dispersion and $v_{\bf k}=\frac{d\epsilon_{\bf k}}{d{\bf k}}$. Note that the energy per particle is $\epsilon_{\bf k}-\mu$. The initial temperature imbalance is characterized by a spatially varying inverse temperature $\beta(x)$ and chemical potential $\mu(x)$.

In the main manuscript we first considered a gas expanding into the vacuum. We model this in $d=1$ by taking $n(x,v,0)=\frac{m}{2\pi}\frac{1}{e^{\beta(mv^{2}/2-\mu)}-1}$ when $x\in A$, and zero outside $A$. For bosons $\mu<0$, so let’s define $\alpha=e^{-\beta\mu}>1$. The particle density at $x=0$ at late time $t\gg 1$ then equals

[TABLE]

Similarly, the energy of bosons $\epsilon_{\bf k}-\mu$ is always positive and nonzero which implies that the late-time behavior of the energy is $E(x=0,t)=\int dv(\epsilon(v)-\mu)n(-vt,0)\sim t^{-1}$. Therefore, whenever a system thermalizes with a powerlaw $t^{-d}$, we will call this ballistic behavior.

Finally, in the main manuscript we show that the correct quantum description of a hot bosonic system $A$ in a cold bosonic bath displays diffusive rather than ballistic behavior. In Fig. B.5 we compare the results of this quantum thermalization to the classical ballistic picture following Eqn. (18) at $\beta_{B}=5$. The classical picture incorrectly yields a ballistic spread, while the exact quantum computation displays diffusive behavior. The diffusive behavior for cold bosonic baths is therefore a true quantum effect.

Appendix C Entropy

Entropy plays an important role in quantum many-body physics. Fortunately, the total entropy of the system, which is of course time-independent, is relatively easy to compute using the modular Hamiltonian,

[TABLE]

For a free system, this implies that the entropy can be expressed in term of the eigenvalues of the Greens function $g_{\alpha}$ as

[TABLE]

where $\eta=\pm 1$ is the sign for bosons/fermions.

To obtain the entanglement entropy of the subsystem $A$, we need the reduced density matrix $\rho_{A}(t)=\mathrm{Tr}\;_{B}\rho(t)$. However, for free systems the entanglement entropy can be computed simply by using Eqn. (23) where $g_{\alpha}$ are now the eigenvalues of the Greens function restricted to subsystem $A$.

Appendix D Thermalization of relativistic fermions

There is an extensive literature on thermal quenches and thermalization of relativistic particles.[1, 2016arXiv160302889C, 3, Lucas:2015uk] In one dimension, the Hamiltonian for relativistic fermions is

[TABLE]

where $\psi_{L,R}(x)$ is the field for left- and right-moving particles, respectively, and $v>0$ is the Fermi velocity. The energy of the right-movers is $\epsilon^{R}_{k}=vk$ and of left-movers is $\epsilon^{L}_{k}=-vk$. The right-moving nature of the $\psi_{R}$ field becomes obvious when one expresses the time evolution of the operator,

[TABLE]

The initial modular Hamiltonian for our hot cloud in $A$ immersed in a bath reads

[TABLE]

Consider only the right-moving particles $h_{R}(x,t)=\psi^{\dagger}_{R}(x,t)iv\partial_{x}\psi_{R}(x,t)$ in subsystem $A$. Under unitary time evolution this segment shifts in its enterity to the right,

[TABLE]

Similarly, the left-movers move to the left under time evolution. This means that after a time $t=L_{A}/v$ there are no remnants of the hot cloud left in the subsystem $A$. The moment the system $A$ has reached a full causal contact with the bath, it is instantaneously thermalized. In Fig. 1 of the main manuscript, middle left, we display the heat profile of such a relativistic system. Notice also that this procedure correctly reproduces the non-equilibrium steady state (NESS) $\mathcal{M}=\beta_{+}H+\beta_{-}P$ as described in Ref. [3].

Appendix E Decay of modular matrix for nonrelativistic particles

In the main manuscript we showed using the continuum approximation $\epsilon_{\bf k}\approx Jk^{2}-\mu$ that the matrix elements of the modular matrix decay ballistically, $\Delta m\sim t^{-d}$. In Fig. E.6 we indeed show that this is the case when computed numerically on large finite size systems. The solid lines are the exact results, and the dashed line is the continuum limit. Note that there are oscillations visible with angular frequency $4J$, which is the result of corrections due to the lattice dispersion $\epsilon_{\bf k}=-2J\sum_{i=1}^{d}\cos k_{i}$. The precise shape and amplitude of these oscillations depend on the precise form of the dispersion.