Evaporative cooling and self-thermalization in an open system of interacting fermions

TL;DR

This paper investigates how an open system of weakly interacting fermions cools through evaporation, demonstrating rapid thermalization to Fermi-Dirac distribution and significant cooling of remaining particles.

Contribution

It introduces a model showing instantaneous thermalization during fermionic evaporation, connecting laboratory cooling techniques with theoretical models including the SYK black hole.

Findings

System quickly reaches thermal equilibrium during evaporation

Occupation numbers follow Fermi-Dirac distribution at all times

Evaporation results in substantial cooling of the system

Abstract

We study depletion dynamics of an open system of weakly interacting fermions with two-body random interactions. In this model fermions are escaping from the high-energy one-particle orbitals, that mimics the evaporation process used in laboratory experiments with neutral atoms to cool them to ultra-low temperatures. It is shown that due to dynamical thermalization the system instantaneously adjusts to the new chemical potential and temperature, so that occupation numbers of the one-particle orbitals always obey the Fermi-Dirac distribution. In this way we are able to describe the evaporation process which leads to a significant cooling of particles remaining inside the system. We also briefly discuss the evaporation process in the SYK black hole model that corresponds to strongly interacting fermions.

Equations

$$n_k=\frac{1}{e^{\beta ({ϵ}_k-\mu )}+1},$$

$${\hat{H}}^{(N)}=\sum _{k=1}^M{ϵ}_k{\hat{c}}_k^{†}{\hat{c}}_k+\epsilon \sum _{ijkl=1}^M{J}_{ij,kl}{\hat{c}}_i^{†}{\hat{c}}_j^{†}{\hat{c}}_k{\hat{c}}_l,$$

$$\frac{d\mathcal{R}}{dt}=-i[\hat{H},\mathcal{R}]+{\mathcal{L}}_{loss}(\mathcal{R}),$$

$${\mathcal{L}}_{loss}(\mathcal{R})=\frac{\gamma }{2}({\hat{c}}_M^{†}{\hat{c}}_M\mathcal{R}-2{\hat{c}}_M\mathcal{R}{\hat{c}}_M^{†}+\mathcal{R}{\hat{c}}_M^{†}{\hat{c}}_M),$$

$$\bar{N}(t)=\sum _{N=1}^{N_0}N{P}_N(t),P_N(t)=\mathrm{Tr}[{\mathcal{R}}^{(N)}(t)]/N,$$

$${\hat{H}}_{eff}^{(N)}=\sum _{k=1}^M\left({ϵ}_k-i\frac{\gamma }{2}{\delta }_{k,M}\right){\hat{c}}_k^{†}{\hat{c}}_k+\epsilon \sum _{ijkl=1}^M{J}_{ij,kl}{\hat{c}}_i^{†}{\hat{c}}_j^{†}{\hat{c}}_k{\hat{c}}_l,$$

$$\bar{N}(t)=N_0\exp \left[-\alpha \bar{N}(t)t\right],$$

$${\hat{H}}_{eff}|{\varphi }_j⟩={\mathcal{E}}_j|{\varphi }_j⟩,{\mathcal{E}}_j=E_j-i\frac{{\Gamma }_j}{2}.$$

$$\bar{N}=\sum _{k=1}^M\frac{1}{e^{\beta ({ϵ}_k-\mu )}+1},\bar{E}=\sum _{k=1}^M\frac{{ϵ}_k}{e^{\beta ({ϵ}_k-\mu )}+1}.$$

Full text

Quantum chaos; semiclassical methods

1 Introduction

The cooling of neutral atoms to micro-Kelvin and further to nano-Kelvin temperatures is one of the most noticeable achievements of the modern physics that opened a new era of quantum technologies [1]. The main method used to cool atoms from micro to nano-Kelvin temperature is the evaporative cooling, where experimental setups are designed to remove the most hot atoms from an atomic cloud. Although being quite successful the method has a drawback that one loses up to 99 percent of atoms to reach the temperature where Bose or Fermi atoms enter the degenerate state.

In this work we extend our previous studies [2] of dynamical (or self-) thermalization in the Two-Body Random Interaction Model (TBRIM) introduced in Refs. [3, 4]. The model describes a system of $N$ spinless fermions distributed over $M$ one-particle orbitals, where fermions have two-body interactions with random matrix elements. We note in passing that at present the limiting case of this model with strongly interacting fermions got a renewed interest in the context of SYK black hole model for theory of quantum gravity [5, 6, 7, 8].

It was shown in Ref. [2, 9] that at moderate interaction strength TBRIM enters the regime of Quantum Chaos111This assumes the interaction strength to exceed some critical value which can be estimated by using Åberg’s criteria [10]. Numerically the transition to Quantum Chaos [11, 12] is observed as a change of the level-spacing statistics from the Poisson distribution to the Wigner-Dyson distribution when the interaction strength exceeds the above critical value [13]. where its eigenstates $|\Psi_{E}\rangle$ become self-thermalized, in spite of the fact that the system is isolated (i.e., not connected to any thermostat). This means, in particular, that occupation numbers $n_{k}=\langle\Psi_{E}|\hat{c}_{k}^{\dagger}\hat{c}_{k}|\Psi_{E}\rangle$ of the orbitals with energies $\epsilon_{k}$ obeys the Fermi-Dirac distribution,

[TABLE]

where the inverse temperature $\beta$ and the chemical potential $\mu$ are uniquely determined by the eigenstate energy $E$ and the number of fermions $N$ in the system. Here we make the system open by introducing absorption of particles from the most upper orbital. It is shown below that TBRIM preserves its self-thermalization property also in the presence of absorption.

Since we remove particles from the high-energy orbital, the open TBRIM mimics the process of evaporative cooling of fermionic atoms. We argue in the work that the understanding of the mechanisms of self-thermalization allows to perform an optimization of the evaporative cooling. It involves, in particular, optimization of the rate at which particles are absorbed in the system. We find that the optimal rate is actually determined by the self-thermalization rate in TBRIM.

2 The model

The Hamiltonian of the closed (isolated) TBRIM reads

[TABLE]

where the orbital energies $\epsilon_{k}$ and the interaction (generally complex) constants $J_{ij,kl}$ are random numbers and we set the dispersion of all random entries to unity. The parameter $\varepsilon$ in the Hamiltonian (2) controls the strength of two-body interactions and the super-index $N$ denotes the number of fermions in the system. The dimension of the Hilbert space is ${\cal N}_{N}=M!/N!(M-N)!$. In what follows we use $M=12$ and $\varepsilon=0.034$. We checked that for this value of $\varepsilon$ the TBRIM is chaotic and self-thermalized for $2\leq N\leq 5$, which are relevant for the numerical simulation reported below.

We describe the evaporation dynamics by solving the master equation on the system density matrix ${\cal R}(t)$,

[TABLE]

[TABLE]

where $\gamma$ is the absorption rate or depletion constant, i.e., the rate at which particles are removed from the most upper one-particle orbital $\epsilon_{M}$. Notice that the density matrix ${\cal R}$ is defined in the extended Hilbert space given by direct sum of $N_{0}+1$ subspaces where $N_{0}=N(t=0)$ is the initial number of fermions. Correspondently, the Hamiltonian $\widehat{H}$ in Eq. (3) has block structure with blocks given by Eq. (2). In what follows we consider $N_{0}=5$ where the total dimension of the Hilbert space is ${\cal N}=1+12+66+220+495+792=3003$. As the initial condition we choose an eigenstate of $\widehat{H}^{(N)}$ with $N=N_{0}$ from the middle of its spectrum. This choice corresponds to infinite effective temperature where occupation numbers of the natural orbitals are approximately equal. The quantities which we calculate are the occupation numbers $n_{k}(t)={\rm Tr}[\hat{c}_{k}^{\dagger}\hat{c}_{k}{\cal R}(t)]$, the mean number of fermions in the system, $\bar{N}(t)=\sum_{k=1}^{M}n_{k}(t)$, and the mean energy, $\bar{E}(t)=\sum_{k=1}^{M}\epsilon_{k}n_{k}(t)$. Typical behavior of $\bar{N}(t)$ and $\bar{E}(t)$ is exemplified in Fig. 1(a,b) where the blue solid and red dashed lines refer to the case of a small $\gamma=0.2$ and a large $\gamma=2$, respectively. In the subsequent two sections we analyze this behavior in some details.

3 Depletion dynamics

We begin with the depletion dynamics. Similar to the total Hamiltonian the density matrix ${\cal R}(t)$ has the block structure where each block is associated with the fixed number of fermions in the system. This gives another expression for the mean number of particles,

[TABLE]

where $P_{N}(t)$ are interpreted as probabilities to find $N$ particles in the system at a given time $t$. Probabilities $P_{N}(t)$ are shown in Fig. 2 in the linear and logarithmic scales for $\gamma=2$. It is seen that $P_{N}(t)$ with $N>1$ show asymptotic exponential decay while $P_{1}(t)$ approaches unity. This is consistent with the expectation that the steady-state solution of the master equation (3) corresponds to a single fermion. It is also seen in Fig. 2 that relaxation to this steady-state is a cascade-like process where ‘children’ cascade takes essentially longer time than ‘parent’ cascade. In the other words, depletion rate decreases proportionally to the number of particles left in the system. In the next paragraph we give a formal explanation for this intuitively expected result and quantify the dependence $\bar{N}(t)$.

One gets a useful insight in the depletion process by analyzing the survival probability, i.e., probability to find the initial number of particles. The survival probability has no parent cascade and can be calculate in a simpler way, namely, by solving the Schrödinger equation with the non-Hermintian Hamiltonian $\widehat{H}_{eff}$,

[TABLE]

which is obtained from the Herminia Hamiltonian (2) by prescribing imaginary energy to the $M$th orbital. Then the norm of the wave function exactly corresponds to probability to find the initial number of particles in the system. The introduced non-Herminia Hamiltonian is related to the problem of quantum chaotic scattering [14, 15, 16, 17, 18]. Considering the Hamiltonian matrix in the Fock basis we find the number of complex diagonal elements to be given by $Q=N{\cal N}_{N}/M$, which can be interpreted as the number of decay channels. It is known that short-time dynamics of the survival probability in chaotic scattering is the exponential decay, $P_{N}(t)=\exp(-\nu t)$, with increment $\nu$ proportional to ratio of the number of channels to the matrix size [17, 18]. Adopting this result to our problem with cascade dynamics (where the number of channels and the matrix size are changing from parent to children cascade) we obtain

[TABLE]

where the coefficient $\alpha$ is a function of the depletion constant $\gamma$. The solution of the nonlinear equation (7) is depicted in Fig. 1(a) by dotted lines. It is in a reasonable agreement with the straightforward simulation of the depletion dynamics.

At the next step we discuss the dependence of the increment $\nu$ for the exponential decay of the survival probability on the depletion constant $\gamma$. One may naively expect that $\nu$ [and, hence, the coefficient $\alpha$ in Eq. (7)] linearly depends on $\gamma$. However, this is true only if the depletion constant is smaller than the rate of self-thermalization, while in the opposite case $\nu$ actually decreases with an increase of $\gamma$. This effect, which is often referred to as the Zeno effect [19], has a simple explanation in terms of the energy spectrum of the Hamiltonian (6),

[TABLE]

Indeed, for small $\gamma$ the complex energies ${\cal E}_{j}$ are located near the real axis, see Fig. 3(a), and the mean value $\bar{\Gamma}$ of the resonance widths is proportional to $\gamma$, $\bar{\Gamma}\sim\gamma$. As $\gamma$ is increased, the cloud of eigenenergies splits into two groups222In the context of quantum chaotic scattering this effect is discussed in Refs. [15, 16]., see Fig. 3(b). For the first group $\bar{\Gamma}$ remains to be proportional to $\gamma$ while for the the second group we have $\bar{\Gamma}\sim 1/\gamma$ , as it is easy to show by using the first-order perturbation theory. Thus, with an increase of the depletion constant $\gamma$ the increment $\nu=\bar{\Gamma}$ shows a maximum at some critical value $\gamma_{cr}$, which for the parameters of Fig. 2 corresponds to $\gamma_{cr}\approx 1.2$. We mention that this critical value of the depletion constant can be used as unambiguous definition of the self-thermalization rate in the considered system of weakly-interacting fermions.

4 Thermalization dynamics

In this section we provide numerically evidence that TBRIM remains to be self-thermalized also in the presence of particle loss. The numerical confirmation of this statement is given below.

First, using numerical data for the mean number of particle $\bar{N}$ and the mean energy $\bar{E}$ depicted in Fig. 1(a,b) we solve the system of two nonlinear algebraic equation on the chemical potential $\mu$ and the inverse temperature $\beta$,

[TABLE]

The solution is shown in the panels (c) and (d) where, as before, the solid and dashed lines refer to $\gamma=0.2$ and $\gamma=2$, respectively. Second, using the obtained $\mu$ and $\beta$ we calculate the occupation numbers according to Eq. (1) and compare them with actual occupation numbers calculated as $n_{k}(t)={\rm Tr}[\hat{c}_{k}^{\dagger}\hat{c}_{k}{\cal R}(t)]$, see Fig. 4. (From now on we focus on the case $\gamma=0.2$.) The observed agreement confirms that we indeed have the true thermalization where the notion of temperature is absolutely meaningful, in spite of the absence of any thermal bath. From Fig.1(d) we clearly see that during evaporation the temperature of the system decreases with time,

In addition to Fig. 4, Fig. 5 shows occupation numbers at the beginning and the end of numerical simulation as the function of orbital energies. Remarkably, occupations of two lowest-energy orbitals get increased. Thus, for a larger system size one may expect that $n_{k}$ of a few lowest-energy orbitals approach unity, i.e., we enter the degenerate state.

5 Case of SYK black hole model

Due to a significant recent interest to the SYK model [5, 6, 7, 8] we briefly discuss the evaporation process in this model which corresponds to vanishing orbital energies $\epsilon_{k}=0$ in the Hamiltonian (2). Clearly, in this case of strongly interacting fermions we cannot appeal to the Fermi-Dirac distribution (which is strictly derived for non-interacting fermions) and, hence, the notion of temperature is not defined. However, we still can address the depletion dynamics.

As for the previous TBRIM case we assume the absorption to take place only from one orbital. Upper panels in Fig. 6 shows probabilities $P_{N}(t)$ for $\gamma=2$ in the linear and logarithmic scales. Similar to the TBRIM case there is an exponential decay of probability to find the number of fermions $N>1$ while the probability $N=1$ approaches unity. The decay rate of the survival probability is seen to be slightly larger than in the TBRIM case, which is consistent with distribution of eigenvalues of the non-Hermitian Hamiltonian – compare Fig. 3(b) and Fig. 6(d). For the sake of completeness we also display in Fig. 6(c) the eigenvalue distribution for $\gamma=0.2$. Similar to the TBRIM case this compact cloud of eigenvalues is found to separate into two clouds at $\gamma\approx 1.2$. Thus for both TBRIM and SYK models the depletion rate is maximized at the same value of $\gamma_{cr}$.

6 Conclusions

We analyzed the process of evaporative cooling in the system of weakly interacting fermions. This process has a competition between depletion, where particles are removed from high-energy one-particle orbitals (generalization of the results onto the case of more than one decaying orbital is straightforward) and self-thermalization, which repopulates these orbitals. We especially stress the importance of the latter process and it is generally a hard problem to find conditions for emergence of self-thermalization in a given system of interacting particles. For the isolated TBRM these conditions have been analyzed in our previous works [2, 9]. In the present work we showed that the self-thermalization in TBRIM takes place also in the presence of particle loss (absorption). Thus the TBRIM becomes a paradigmatic model for investigation of different aspects of self-thermalization in closed and open systems.

We studied the depletion dynamics in the open TBRM by mapping the problem to quantum chaotic scattering. It was shown, in particular, that the number of lost particles is not a monotonic function of the depletion constant $\gamma$ but has a pronounced maximum at some $\gamma_{cr}$. Since the depletion constant can be varied in laboratory experiments, this result suggests a method for measuring the self-thermalization rate by finding the critical $\gamma_{cr}$ where the increment $\nu$ for exponential decay of the survival probability is maximized.

Since we remove particles from the high-energy orbitals, the depletion process results in a decrease of the system energy and, as a consequence, in the temperature drop for the remaining particles. We stress that efficiency of this cooling mechanism is not just proportional to particle loss, as one might naively expect, but is an involved function of the depletion constant $\gamma$ (which can be varied in a laboratory experiment), the self-thermalization rate $\gamma_{cr}$ (which is an internal property of the system), and the duration of the evaporation process (which can be limited by some reason). An example is given in Fig. 1(a,d). It is seen that within the same time interval we reached lower temperature for $\gamma=0.2$ than for $\gamma=2$, in spite of the fact that in the former case we lost less particles than in the later case.

7 Acknowledgements

For DLS this work was supported in part by the Programme Investissements d’Avenir ANR-11-IDEX-0002-02, reference ANR-10-LABX-0037-NEXT (project THETRACOM).