Simulating quantum thermodynamics of a finite system and bath with variable temperature

TL;DR

This paper models a finite quantum bath with variable temperature to simulate heat exchange and thermalization in small quantum systems, revealing finite-size effects and deviations from classical thermodynamic relations.

Contribution

It introduces a finite, variable-temperature quantum bath with non-identical oscillators for simulating quantum thermodynamics, highlighting finite-size effects and deviations from classical behavior.

Findings

Finite baths exhibit higher temperatures than infinite baths at the same energy.

The temperature approaches the infinite bath limit as the number of oscillators increases.

Significant deviations from classical energy-temperature relations are observed in small quantum systems.

Abstract

We construct a finite bath with variable temperature for quantum thermodynamic simulations in which heat flows between a system $\mathcal{S}$ and the bath environment $\mathcal{E}$ in time evolution of an initial $\mathcal{SE}$ pure state. The bath consists of harmonic oscillators that are not necessarily identical. Baths of various numbers of oscillators are considered; a bath with five oscillators is used in the simulations. The bath has a temperature-like level distribution. This leads to definition of a system-environment microcanonical temperature $T_\mathcal{SE}(t)$ which varies with time. The quantum state evolves toward an equilibrium state which is thermal-like, but there is significant deviation from the ordinary energy-temperature relation that holds for an infinite quantum bath, e.g. an infinite system of identical oscillators. There are also deviations from the Einstein quantum heat capacity. The temperature of the finite bath is systematically greater for a given energy than the infinite bath temperature, and asymptotically approaches the latter as the number of oscillators increases. It is suggested that realizations of these finite-size effects may be attained in computational and experimental dynamics of small molecules.

Tables

Table 1. Table 1: Oscillator frequencies in the five harmonic oscillator environment shown to six decimal places.

$\hslash {\omega }_1$	$\hslash {\omega }_2$	$\hslash {\omega }_3$	$\hslash {\omega }_4$	$\hslash {\omega }_5$
0.620 246	0.735 401	1.146 315	1.316 886	1.453 415

Table 2. Table 2: Energy and temperature data for Fig. 9 . The energies E = ⟨ H ^ ⟩ 𝐸 delimited-⟨⟩ ^ 𝐻 E=\langle\hat{H}\rangle are from the full Hamiltonian in Eq. 1 and the T S ​ E ​ ( E ) subscript 𝑇 𝑆 𝐸 𝐸 T_{SE}(E) were calculated from Eq. 23 . The average bath-oscillator energies ⟨ E o ​ s ​ c ⟩ = E ℰ / η delimited-⟨⟩ subscript 𝐸 𝑜 𝑠 𝑐 subscript 𝐸 ℰ 𝜂 \langle E_{osc}\rangle=E_{\mathcal{E}}/\eta were averaged over the same time window 30 < t ≤ 60 30 𝑡 60 30<t\leq 60 as the system probabilities in Fig. 9 and the infinite bath T 𝑇 T were calculated from Eq. 14 with ⟨ n o ​ s ​ c ⟩ = ⟨ E o ​ s ​ c ⟩ delimited-⟨⟩ subscript 𝑛 𝑜 𝑠 𝑐 delimited-⟨⟩ subscript 𝐸 𝑜 𝑠 𝑐 \langle n_{osc}\rangle=\langle E_{osc}\rangle .

State	$E$	$T_{𝒮ℰ}$	$⟨E_{osc}⟩$	$T$ (Eq. 14)
(a)	4.148	1.422	0.750 $\pm$ 0.005	1.180 $\pm$ 0.006
(b)	6.118	1.913	1.121 $\pm$ 0.003	1.568 $\pm$ 0.003
(c)	8.099	2.406	1.499 $\pm$ 0.002	1.957 $\pm$ 0.002

Equations

$$\hat{H}={\hat{H}}_{\mathcal{S}}+{\hat{H}}_{\mathcal{E}}+{\hat{H}}_{\mathcal{S}\mathcal{E}}$$

$$0\le E_n+E_{ϵ}\le 13,$$

$$⟨n|{\hat{H}}_{\mathcal{S}}|n⟩=\hslash {\omega }_{\mathcal{S}}n,$$

$${\rho }_{Ein}(\eta ,n_{tot})=\frac{(\eta -1+n_{tot})!}{(\eta -1)!n_{tot}!},$$

$${\hat{H}}_{\mathcal{E}}=\sum _{osc=1}^{\eta }{\hat{H}}_{osc},$$

$$⟨n_{osc}|{\hat{H}}_{osc}|n_{osc}⟩=\hslash {\omega }_{osc}{n}_{osc},$$

$${\rho }_{\mathcal{E}}(E_{\mathcal{E}})=\frac{\Gamma (\eta +E_{\mathcal{E}})}{\Gamma (\eta )\Gamma (E_{\mathcal{E}}+1)},$$

$$\sqrt[\eta ]{\prod _{osc=1}^{\eta }\hslash {\omega }_{osc}}=\hslash \omega =1,$$

$${\rho }_{{\hat{H}}_{\mathcal{E}}}\approx {\rho }_{\mathcal{E}}={\rho }_{Ein}$$

$${\rho }_{{\hat{H}}_{\mathcal{E}}}\approx {\rho }_{\mathcal{E}}$$

$$\frac{1}{T}=\frac{\partial S}{\partial E}$$

$$\frac{1}{T_{\mathcal{E}}}=\frac{d{\rho }_{\mathcal{E}}/d{E}_{\mathcal{E}}}{{\rho }_{\mathcal{E}}}.$$

$$\frac{1}{T_{\mathcal{E}}}=\psi (E_{\mathcal{E}}+\eta )-\psi (E_{\mathcal{E}}+1)=\sum _{m=1}^{\eta -1}\frac{1}{E_{\mathcal{E}}+m},$$

$$⟨n_{osc}⟩=\frac{1}{e^{1/T}-1}$$

$$\underset{⟨n_{osc}⟩\to \infty }{\lim }T=⟨n_{osc}⟩+\frac{1}{2}=⟨E_{osc}^{(+zp)}⟩,$$

$$\underset{E_{\mathcal{E}}\to \infty }{\lim }{T}_{\mathcal{E}}=\frac{E_{\mathcal{E}}+\eta /2}{\eta -1}=\left(⟨E_{osc}⟩+\frac{1}{2}\right)\frac{\eta }{\eta -1},$$

$$C=\frac{d⟨E_{osc}⟩}{dT}$$

$$\frac{1}{T_{\mathcal{S}\mathcal{E}}(E)}=\frac{d{\rho }_{\mathcal{S}\mathcal{E}}/dE}{{\rho }_{\mathcal{S}\mathcal{E}}}.$$

$${\rho }_{\mathcal{S}\mathcal{E}}(E)=\sum _n{\rho }_{\mathcal{E}}(E-E_n).$$

$$\frac{1}{T_{\mathcal{S}\mathcal{E}}(E)}=\sum _n\frac{d{\rho }_{\mathcal{E}}(E-E_n)/dE}{{\sum }_m{\rho }_{\mathcal{E}}(E-E_m)}.$$

$$\frac{1}{T_{\mathcal{S}\mathcal{E}}(E)}=\sum _n\frac{{\rho }_{\mathcal{E}}(E-E_n)}{{\sum }_m{\rho }_{\mathcal{E}}(E-E_m)}\frac{1}{T_{\mathcal{E}}(E-E_n)}.$$

$$\frac{{\rho }_{\mathcal{E}}(E-E_n)}{{\sum }_m{\rho }_{\mathcal{E}}(E-E_m)}=p_{micro}(E_n).$$

$$\frac{1}{T_{\mathcal{S}\mathcal{E}}(E)}=\sum _n\frac{p_{micro}(E_n)}{T_{\mathcal{E}}(E-E_n)}={⟨\frac{1}{T_{\mathcal{E}}(E-E_n)}⟩}_{micro}$$

$$\frac{1}{T_{\mathcal{S}\mathcal{E}}(E,t)}=\sum _n\frac{{\rho }_{\mathcal{S}}^{n,n}(t)}{T_{\mathcal{E}}(E-E_n)}={⟨\frac{1}{T_{\mathcal{E}}(E-E_n)}⟩}_{{\hat{\rho }}_{\mathcal{S}}(t)}$$

$$|{\Psi }_{n_0}⟩=|n_0⟩|{ϵ}_0⟩,$$

$$|{ϵ}_0⟩\sim \sum _{ϵ}\exp \left(-\frac{(E_{ϵ}-E_{{ϵ}_0}{)}^2}{2{\sigma }_{\mathcal{E}}^2}\right)|ϵ⟩,$$

$$E_{{ϵ}_0}=E_0-E_{n_0}$$

$$⟨n|⟨ϵ|\hat{R}|{ϵ}^{\prime }⟩|n^{\prime }⟩=R_{nϵ,n^{\prime }{ϵ}^{\prime }}.$$

$$p(X_{nϵ,n^{\prime }{ϵ}^{\prime }})\sim e^{-X_{nϵ,n^{\prime }{ϵ}^{\prime }}^2/2{\sigma }^2},$$

$$⟨n|⟨ϵ|{\hat{H}}_{\mathcal{S}\mathcal{E}}|{ϵ}^{\prime }⟩|n^{\prime }⟩=k{R}_{nϵ,n^{\prime }{ϵ}^{\prime }}\exp \left(-{\gamma }_{\mathcal{S}}|\Delta n|-{\gamma }_{\mathcal{E}}\sum _{osc=1}^{\eta }|\Delta n_{osc}|\right)$$

Full text

I introduction

This paper considers computational simulation of a process of energy flow as a quantum system becomes entangled with a very small temperature bath. In the corresponding “classical” thermodynamic system, we would have an idea of a variable temperature as energy flows into the finite bath. Here we ask, does a simulacrum of thermodynamic behavior emerge when we make the bath very small? Do reasonable ideas of a variable temperature hold, and is there something akin to thermal equilibrium with a Boltzmann distribution? We will find that with a very small “thermal” environment, as small as five oscillators, it is possible to get behavior that is very much like thermodynamic behavior. On the other hand, anomalies are observed related to the notion of temperature with the small bath. The work here builds on earlier simulations with a cruder, constant temperature bath Barnes and Kellman (2013); Barnes et al. (2018); Lotshaw and Kellman (2019); Borowski et al. (2003); Silvestri et al. (2014); Esposito and Gaspard (2003). Questions of variable temperature in a very small quantum thermodynamic system and bath are of more than abstract interest. Our simulations may not be too much simpler than what is called for in problems of practical import. Quantum nanodevices can be imagined whose performance may depend on considerations similar to those here. Similar in spirit to the approach taken here, quantum thermalization behavior of a pure quantum state has recently been observed experimentally in Bose-Einstein condensates containing as few as six-atoms Kaufman et al. (2016). Recently Pérez and Arce (2018); Leitner (2015, 2018), work on molecular “quantum chaos” is being conceptualized as a venue for the exploration of contemporary ideas about the foundations of quantum thermodynamics, to which we turn next.

There have been a variety of simulations of quantum thermodynamic processes, including the very basic elementary process of heat flow into a bath Barnes and Kellman (2013); Barnes et al. (2018); Lotshaw and Kellman (2019); Borowski et al. (2003); Silvestri et al. (2014); Esposito and Gaspard (2003). These have been successful in recovering standard thermodynamic behavior, with attainment of thermal equilibrium and a Boltzmann distribution for the system, with a properly behaving temperature. However, these investigations have used rather simple models of the temperature bath, sometimes with a grossly discrete model of energy levels Barnes and Kellman (2013); Barnes et al. (2018); Borowski et al. (2003), in others with an approximation to continuous levels in the bath Lotshaw and Kellman (2019); Silvestri et al. (2014); Esposito and Gaspard (2003), but always to our knowledge with a model of an effectively infinite bath with fixed temperature in mind. Usually also, a very simple coupling between system and environment is assumed, typically, a random matrix coupling without significant structure. Paralleling (and sometimes preceding) these simulations, there has been a great deal of work Barnes et al. (2018); Lotshaw and Kellman (2019); Leitner (2015, 2018); Rigol et al. (2008); Deutsch (2018, 1991); D’Alessio et al. (2016); Tasaki (1998); Gemmer et al. (2009); Popescu et al. (2006); Linden et al. (2009); Goldstein et al. (2006); Goldstein et al. (2010a, b); Goldstein et al. (2015); von Neumann (2010); Reimann (2008, 2016); Esposito et al. (2010); Polkovnikov (2011); Han and Wu (2015); Kak (2007); Reeb and Wolf (2014); Xu et al. (2014); Logan and Wolynes (1990) examining theoretical foundations of quantum thermodynamics. Generally, this has focused on the large $N$ limit of quantum entangled systems. In our simulations here the focus is rather on the extent to which thermodynamic-like behavior persists as the total system becomes very small. There have been simulations examining ergodicity and energy flow in small total systems Leitner (2015, 2018); Rigol et al. (2008); Deutsch (2018); Bigwood and Gruebele (1995); Gruebele (2003), but these have not involved the type of variable temperature analysis that is our focus here. We construct a finite, variable temperature bath, also making use of a structured coupling which is far more selective than the random matrix coupling used in many earlier simulations. We will find that we can build a simulation model with features very much like a variable temperature and thermalization, but with significant anomalies due to the finite bath, with some challenges to overcome having to do with the nature of the coupling.

As noted briefly above, and in more detail in the concluding section, there are real molecular systems that could be considered as laboratories for “post-classical” thermodynamic effects. Consideration of small size is a recent “dimension” of quantum thermodynamics beyond that introduced long ago with the advent of quantum levels. A third innovation might come with novel effects from combining quantum time evolution with multiple small baths of the kind developed here for a single bath.

II Model System-Environment “Universe”

In this section, we detail the system and environment in our model; we treat the system-environment interaction separately, in Sections V and VI .

We will deal with a total system or “universe” pure state for a coupled and entangled system and environment, or temperature bath. The total Hamiltonian includes system $\mathcal{S}$, environment $\mathcal{E}$, and interaction \mathcal{S}$$\mathcal{E} components

[TABLE]

For the basis set we will use a truncation of the full \mathcal{S}$$\mathcal{E} tensor product basis to a subset that contains all of the \mathcal{S}$$\mathcal{E} basis states $|n\rangle\otimes|\epsilon\rangle$ in the energy range

[TABLE]

similar to the “thermal basis” described in Ref. Lotshaw and Kellman (2019). The numerical convergence with this basis will be discussed in Section VI. Time evolution of the pure \mathcal{S}$$\mathcal{E} state $|\Psi\rangle$ is carried out by numerically diagonalizing $\hat{H}$ and then calculating a series of timesteps using the Schrödinger equation $|\Psi(t)\rangle=\exp(-i\hat{H}t)|\Psi(0)\rangle$ ($\hbar=1)$. In this section we will develop the system and environment basis sets and Hamiltonians $\hat{H}_{\mathcal{S}}$ and $\hat{H}_{\mathcal{E}}$; later sections develop $\hat{H}_{\mathcal{SE}}$.

The system Hamiltonian consists of a set of five evenly spaced levels

[TABLE]

with frequency $\omega_{\mathcal{S}}=$ 0.5 and quantum number $n=0,1,...,4$. These choices of $\omega_{\mathcal{S}}$ and $n$ give a maximum system energy $E_{\mathcal{S}}^{max}=2$ that is reasonably small compared to the initial \mathcal{S}$$\mathcal{E} state total energies we will consider in this paper $\langle\hat{H}\rangle\gtrsim 4$, where $\hat{H}$ is the total Hamiltonian of Eq. 1. With larger $E_{\mathcal{S}}^{max}$ we have found that it is more difficult to get good system thermalization, since very few environment levels are paired with the highest energy system levels at the total energy $\langle\hat{H}\rangle$ when $E_{\mathcal{S}}^{max}\approx\langle\hat{H}\rangle$. This choice of $\omega_{\mathcal{S}}$ and $n=0,1,...,4$ ensures that there is always a fair amount of energy in the environment, so that it can act properly as a heat bath to the system in our simulations.

We want to have an environment or bath $\mathcal{E}$ with certain properties more general than in earlier work Barnes and Kellman (2013); Barnes et al. (2018); Lotshaw and Kellman (2019); Borowski et al. (2003); Silvestri et al. (2014); Esposito and Gaspard (2003), and more similar to real physical systems. We want the temperature to vary with energy, instead of being fixed. We would also like for the energy and temperature to be close to proportional, $T\sim E$, to the extent possible in a finite model, and exactly so in the limit of a large bath. Furthermore, we may want the bath to have some significant structure, so that the couplings might also have some structure, unlike the abstract undefined environment levels with random couplings used earlier. To do all of these things, we will construct the bath as a collection of oscillators.

Consider first a set of degenerate oscillators with equal frequencies and level spacings $\ \hbar\omega=1$. This “Einstein heat capacity” system has the well known degeneracy pattern and density of states

[TABLE]

where $\rho_{Ein}(\eta,n_{tot})$ is the number of ways to distribute $n_{tot}$ total energy quanta into $\eta$ oscillators. A more physically realistic model will generalize to oscillators of different frequencies, so as to obtain something resembling a continuous distribution of levels, while approximately maintaining the overall pattern of Eq. 4. To this end, we will extend the distribution $\rho_{Ein}$ to variable frequencies and energies using a continuous function $\rho_{\mathcal{E}}$ that interpolates between the discrete points in Eq. 4. Then, we will devise a set of distinct harmonic oscillator frequencies $\{\omega_{osc}\}$ that approximates the continuous distribution. The total environment Hamiltonian is expressed as the sum of oscillator Hamiltonians

[TABLE]

where the $\hat{H}_{osc}$ have energy eigenvalues

[TABLE]

where $n_{osc}$ is the quantum number of a given oscillator. We will analyze the density of states $\rho_{\hat{H}_{\mathcal{E}}}$ of the Hamiltonian $\hat{H}_{\mathcal{E}}$, finding good agreement with the continuous density $\rho_{\mathcal{E}}$, and then analyze the temperature dependence of the model.

We begin by developing a continuous density function $\rho_{\mathcal{E}}$ in place of the highly degenerate density of Eq. 4. The most straightforward way to do this is to replace the factorials in (4) with Gamma functions

[TABLE]

where the discrete number of total quanta $n_{tot}$ has been replaced by a continuous environment energy $E_{\mathcal{E}}$. The $\Gamma$ function extends the density to non-integer values of the energy $E_{\mathcal{E}}$, and agrees with the original density $\rho_{Ein}$ at integer $E_{\mathcal{E}}$ $=n_{tot}$, since for example $\Gamma(E_{\mathcal{E}}+1)=E_{\mathcal{E}}!$ $=n_{tot}!$ when $E_{\mathcal{E}}$ $=n_{tot}$ is an integer. The top of Fig. 1 shows how the continuous density $\rho_{\mathcal{E}}$ extends the degenerate oscillator density $\rho_{Ein}$ to non-integer $E_{\mathcal{E}}$.

The next step is to devise a set of oscillator frequencies for the Hamiltonian $\hat{H}_{\mathcal{E}}$ in Eq. 5 with a density $\rho_{\hat{H}_{\mathcal{E}}}$ that follows the interpolating function $\rho_{\mathcal{E}}$. An $\eta=5$ oscillator bath will be used for the simulations. This value of $\eta$ is large enough to give a density of states with an exponential-like dependence on energy, which will be imperative for Boltzmann thermalization of the system $\mathcal{S}$, but also small enough to make the computations tractable. The frequencies are generated as random numbers, to make the bath generic. We first tried generating random numbers $0.5\leq\hbar\omega_{osc}\leq 1.5$ then rescaling the  $\hbar\omega_{osc}$ so that their average was the same as the degenerate oscillator frequency $\hbar\omega=1$ seen in the top of Fig. 1. However, when constructing the Hamiltonian $\hat{H}_{\mathcal{E}}$ in Eq. 5 using these frequencies, it was found that the resulting density of states $\rho_{\hat{H}_{\mathcal{E}}}$ was always greater than the desired $\rho_{\mathcal{E}}$ of Eq. 7. Instead, good agreement $\rho_{\hat{H}_{\mathcal{E}}}\approx\rho_{\mathcal{E}}$ is consistently found by rescaling the random  $\hbar\omega_{osc}$ values according to their geometric mean,

[TABLE]

as discussed in detail shortly. Eq. 8 sets the unit of energy in this paper and also sets the relationship between the collection of variable frequencies $\{\hbar\omega_{osc}\}$ and the degenerate oscillator frequency $\hbar\omega$ assumed in connection with Eq. 4. The relation Eq. 8 has previously been noted by Landau and Lifshitz Landau and Lifshitz (1980) where it was also found to give the necessary link between variable and fixed frequency oscillators in a different context.

The $\hat{H}_{\mathcal{E}}$ that we use with Eq. 5 uses the frequencies given in Table 1 that come from randomly chosen values that have been rescaled according to Eq. 8. The results are robust for other choices of random and rescaled $\{\hbar\omega_{osc}\}$. The density of states $\rho_{\hat{H}_{\mathcal{E}}}$ for this set of frequencies is shown in the histogram boxes in the bottom of Fig. 1, and is in excellent agreement with $\rho_{\mathcal{E}}$ of Eq. 7. Recall that $\rho_{\mathcal{E}}$ also agrees with the fixed frequency $\rho_{Ein}$ as seen in the top of Fig. 1. This demonstrates that Eq. 8 gives the desired correspondence between the densities of states for the variable and identical frequency oscillators:

[TABLE]

at integer energies $E_{\mathcal{E}}=n_{tot}$ and

[TABLE]

at non-integer energies (where the single-frequency $\rho_{Ein}$ is undefined in Eq. 4). The correspondence between the somewhat random $\rho_{\hat{H}_{\mathcal{E}}}$ and the well-controlled, analytical $\rho_{\mathcal{E}}$ will allow us to determine analytical temperature relationships for our oscillator bath using the relatively simple function $\rho_{\mathcal{E}}$. This is developed in the next section.

III Temperature

This rather involved section addresses key questions about the “thermal” character introduced by the small finite bath in our model. Does the standard infinite bath relation $E\sim T$ hold at high energy? What is the low temperature behavior of the finite bath? While sensible notions of temperature will emerge, we will also see that there are anomalies in both of these aspects, related to the finite size of the bath.

We usually think of temperature in terms of a microcanonical ensemble with a very large, effectively infinite bath, so that the temperature is constant. The temperature comes from the standard relation

[TABLE]

applied to the total system+environment ${\mathcal{SE}}$ as the density of states is varied with energy. In the situation envisaged in Fig. 2, we start by thinking instead of a temperature $T_{\mathcal{E}}$ for the bath environment initially in isolation from the system. There are a multiplicity of initial separate system-bath combinations, each with the same total energy $E$; an example is the red \mathcal{S}$$\mathcal{E} state pair in the left of Fig. 2. Each ${\mathcal{SE}}$ combination has its own initial system energy $E_{\mathcal{S}}$, bath energy $E_{\mathcal{E}}$, and bath temperature $T_{\mathcal{E}}$. The bath temperature $T_{\mathcal{E}}$ is based on a fixed $E_{\mathcal{E}}$ microcanonical energy that is defined only before the interaction with the system has begun – the system in our simulations starts in a single zero-order state – so there is no meaningful independent system temperature. Then, heat flows between system and bath, leading to a finite change in a temperature that we want to be defined for the final equilibrium state, and perhaps in between as well. The final temperature $T_{\mathcal{SE}}$ after the heat flow comes from the microcanonical ensemble for the total system \mathcal{S}$$\mathcal{E}, which consists of the union of all the system-bath sub-ensembles, all with total \mathcal{S}$$\mathcal{E} energy $E$, as in the right of Fig. 2. An interesting relation Eq. 23 will be found to hold between the inverse temperature $1/T_{\mathcal{SE}}$ of the complete ensemble of the \mathcal{S}$$\mathcal{E} total system, and the average of the inverse temperatures $1/T_{\mathcal{E}}$ of the baths of the sub-ensembles. In fact, it will be possible to define a time-varying “master temperature” $T_{\mathcal{SE}}(t)$ in Eq. 24 for the time-dependent intermediate state $|\Psi(t)\rangle$ in the equilibration process. Thus, we will obtain a satisfying unified description of all the possible processes of the type in Fig. 2.

III.1 Temperature for Initial Isolated Environment

First, we develop the temperature $T_{\mathcal{E}}$ for a finite environment that is thermally isolated from the system. (This will turn out to be the initial state temperature in the time-dependent temperature $T_{\mathcal{SE}}(t)$ to be developed in Section III.3.) We will compare this finite bath to an infinite “true” temperature bath of infinitely many oscillators. The system is in a single zero-order initial state $n_{0}$, corresponding to our initial state in Fig. 2. The total energy is $E$, the system has energy $E_{\mathcal{S}}=E_{n_{0}}$, and the environment has energy $E_{\mathcal{E}}=E-E_{\mathcal{S}}$. The temperature is defined using the standard thermodynamic relation of Eq. 11. This is evaluated using the Boltzmann entropy $S=k_{B}\ln W(n_{0},E)$, with $W(n_{0},E)$ the number of \mathcal{S}$$\mathcal{E} states $|n_{0},\epsilon\rangle$ in a microcanonical energy shell $[E-\delta E/2,E+\delta E/2]$, again with the system in the level $n_{0}$. Since $n_{0}$ is fixed, $W(E)=\rho_{\mathcal{E}}(E_{\mathcal{E}})\delta E$ is just the number of environment states, where $\rho_{\mathcal{E}}$ in Eq. 7 is the smoothed continuous density function describing the density of discrete states in our Hamiltonian $\rho_{\hat{H}_{\mathcal{E}}}$, following Eqs. 9 and 10. The initial temperature is then related only to the environment, and we will label it $T_{\mathcal{E}}$, and rewrite it in terms of the density $\rho_{\mathcal{E}}$ as

[TABLE]

Using Eq. 7 for $\rho_{\mathcal{E}}$ then gives

[TABLE]

where $\psi(x)=(d\Gamma(x)/dx)/\Gamma(x)$ is the digamma function. The last equality comes analytically from $\eta-1$ applications of the recurrence relation wol $\psi(x)=\psi(x-1)+1/(x-1)$ to the term $\psi(E_{\mathcal{E}}+\eta)$.

It is not clear just from looking at Eq. 13 how our temperature $T_{\mathcal{E}}$ for the finite bath will behave in comparison to standard temperature-energy relations involving an infinite fixed-temperature bath. In the next two subsections we will make this comparison, using the paradigmatic standard of an average oscillator in an infinite oscillator bath. Section III.1.1 will discuss the convergence of $T_{\mathcal{E}}$ from Eq. 13 to the standard temperature-energy relation as the size of the bath is increased, with convergence to the high energy relation $T\sim E$. Section III.1.2 will discuss deviations related to the finite size of the bath, including deviations from $T=0$ at low energy, and deviations in the heat capacity even at high energy.

III.1.1 Comparison of finite and infinite bath: energy-temperature relation

The heat bath described above is a finite collection of oscillators. We will compare this to a true temperature bath consisting of an infinite collection of oscillators. For this, we use the energy-temperature relation from Einstein and Planck for a harmonic oscillator in an infinite temperature bath:

[TABLE]

($\hbar\omega=1$ and $k_{B}=1$) , where $\langle n_{osc}\rangle$ is the expected number of energy quanta in the oscillator. (This relation was obtained by Einstein in his heat capacity model Einstein (1989) by treating a solid as a collection of identical oscillators in an exterior temperature bath using the canonical ensemble. The result is the same regardless of the ensemble setup, microcanonical or canonical.) We will find that our $T_{\mathcal{E}}$ for the finite bath behaves much like a standard temperature, but also has significant differences from the Einstein relation Eq. 14, leading also to deviations in the heat capacity from the Einstein model. However, we also find that $T_{\mathcal{E}}$ agrees properly with Eq. 14 in the limit of a large number of oscillators. The development is based on the correspondence $\rho_{\mathcal{E}}\approx\rho_{\hat{H}_{\mathcal{E}}}$ in Eqs. 9 and 10, recalling the remarks there about the analytical function $\rho_{\mathcal{E}}$

These relationships are represented in Fig. 3 and later for the heat capacity in Fig. 4. It will be instructive to consider the total energy of the “Einstein oscillator” including both energy quanta and the zero-point energy, $\langle E_{osc}^{(+zp)}\rangle=\langle n_{osc}\rangle+1/2$. The blue curve in Fig. 3 shows the relationship between $\langle E_{osc}^{(+zp)}\rangle$ and temperature based on Eq. 14. The curve begins at the zero-point energy at $T=0$, then quickly approaches the well-known quantum equipartition relation

[TABLE]

shown by the green line in the background of the figure.

For comparison, Fig. 3 also shows the relationship between $\langle E_{osc}\rangle+1/2$ and $T_{\mathcal{E}}$ for finite oscillator baths with various $\eta$, again, based on the correspondence $\rho_{\mathcal{E}}\approx\rho_{\hat{H}_{\mathcal{E}}}$ in Eqs. 9 and 10. The average energy per oscillator from energy quanta $\langle E_{osc}\rangle\equiv E_{\mathcal{E}}/\eta$ is the analog for our bath of $\langle n_{osc}\rangle$ for the Einstein oscillator in Eqs. 14 and 15. The quantity 1/2 then shifts this up by the Einstein oscillator zero-point energy to allow for a direct comparison in the figure between our $T_{\mathcal{E}}$ and the temperature in the Einstein model. In general, the exact zero-point energy in our model will not be 1/2 in our units (unlike the Einstein model), but will instead depend on the frequencies of the oscillators. Here, the $1/2$ is an arbitrary added quantity for the finite baths, inserted for comparison to the Einstein bath.

For the $\eta=5$ bath we use for our simulations, shown by the black solid curve, the temperature behavior is significantly different than the blue infinite bath curve. As we increase the number of oscillators $\eta$ we find that the curves get closer to the standard blue curve for an infinite bath. For example, the dashed-double-dotted red line for $\eta=500$ oscillators rests on top of the blue line for the infinite bath $T$. The convergence towards Eq. 14 with increasing $\eta$ confirms that our temperature gives the standard relation for an infinite bath in the thermodynamic limit $\eta\to\infty$, as expected with a reasonable temperature definition. With this in mind, we next discuss in more detail the much more interesting question of anomalies in temperature behavior associated with small number of oscillators $\eta$ in the finite bath.

III.1.2 Anomalous temperature behavior associated with a very small bath

The very small size of the $\eta=5$ bath leads to anomalous temperature behavior at both high and low energies, as seen in Fig. 3. First, consider the behavior of $T_{\mathcal{E}}$ at low energies. Recall that we treat this as a continuous variable that will be related to the continuous variable $E_{E}$ in Eq. 13. The temperatures for all of the finite $\eta$ oscillator baths in Fig. 3 are nonzero at the minimum value of energy 1/2 in the figure (when $E_{\mathcal{E}}=0$ in Eq. 13, the rationale for the 1/2 being that given in the last subsection). The non-zero minimum temperatures seem to be an unavoidable consequence of combining a finite bath with the standard temperature definition Eq. 12. The temperature is only zero when $d\rho_{\mathcal{E}}/dE_{\mathcal{E}}=\infty$ in Eq. 12 – an evidently impossible condition for a finite bath with a limited number of states. However, as seen in Fig. 3, the curves for increasing $\eta$ converge to the standard infinite bath relation in which $T=0$ at the minimum energy 1/2.

At high energy, $T_{\mathcal{E}}$ approaches the asymptotic relation

[TABLE]

where again $\langle E_{osc}\rangle=E_{\mathcal{E}}/\eta$ refers to the average energy per non-identical oscillator of the finite bath, although it also applies to an infinite “Einstein bath” of identical oscillators. Eq. 16 comes from the analytical limit of the right-hand side of Eq. 13, which we evaluated using Mathematica. Eq. 16 differs from the high-energy Einstein relation Eq. 15 by the factor of $\eta/(\eta-1)$. This difference is negligible in the thermodynamic limit $\eta\to\infty$ but very significant for small $\eta$, as seen by the differing slopes for the solid black and blue lines in Fig. 3 at high energy.

The differing slopes correspond to a difference in heat capacities

[TABLE]

between the different temperature-energy relations. The heat capacities for all of the temperature-energy curves in Fig. 3 are plotted in Fig. 4. The heat capacity curves are similar to the standard Einstein behavior at low temperature, but they are systematically lower at high temperature, where they approach asymptotic values $C\to(\eta-1)/\eta<1$, less than both the Einstein relation and the standard equipartition result.

We will find in Section VII that the anomalous temperature behavior seen in Fig. 3 is critical in obtaining the correct thermalized Boltzmann distribution for the system: the anomalous scaling behavior $\sim\eta/(\eta-1)$ in the figure must be taken into account to correctly describe the equilibrium $\mathcal{S}$ Boltzmann distribution and the \mathcal{S}$$\mathcal{E} thermodynamic behavior.

III.2 System-Environment Microcanonical Temperature

We now consider the equilibrium \mathcal{S}$$\mathcal{E} state and the temperature $T_{\mathcal{SE}}$ for the complex entangled state $|\Psi(t)\rangle$ shown schematically in the right of Fig. 2; this will be the equilibrium value of the time-dependent temperature $T_{\mathcal{SE}}(t)$ to be developed in Section III.3.

$T_{\mathcal{SE}}$ is defined following the same reasoning leading to Eq. 12, giving

[TABLE]

To evaluate the temperature we will examine $\rho_{\mathcal{SE}}$ as the density of zero-order states, just as we did for the isolated bath temperature $\rho_{\mathcal{E}}$. While there is some arbitrariness in doing this now with $\rho_{\mathcal{SE}}$, it is operationally simple, and seems at least as reasonable a choice as other possibilities. It is consonant with what we have done with $\rho_{\mathcal{E}}$, and will lead to the simple result Eq. 23.

The total density of \mathcal{S}$$\mathcal{E} zero-order states at energy $E$ has contributions from all of the \mathcal{S}$$\mathcal{E} state pairs that are in the microcanonical energy shell $E-\delta E/2\leq E_{\mathcal{S}}+E_{\mathcal{E}}\leq E+\delta E/2$, that is, each of the \mathcal{S}$$\mathcal{E} state pairs shown schematically in Fig. 2. The total density of \mathcal{S}$$\mathcal{E} states is the sum of bath densities that pair with each system level $n$ at the total energy $E=E_{\mathcal{E}}+E_{n}$,

[TABLE]

The \mathcal{S}$$\mathcal{E} temperature can then be written as

[TABLE]

The derivatives can be rewritten in terms of $\rho_{\mathcal{E}}$ and $T_{\mathcal{E}}$ using Eq. 12, giving

[TABLE]

The fraction involving the densities gives the number of microcanonical states with the system in the level $E_{n}$ relative to the total number of microcanonical states. This is simply the microcanonical probability of the system level $E_{n}$,

[TABLE]

Putting this into Eq. 21 gives the simple result

[TABLE]

Eq. 23 says that the reciprocal temperature $1/T_{\mathcal{SE}}$ for the full \mathcal{S}$$\mathcal{E} microcanonical ensemble is simply the average of the reciprocal environment temperatures $1/T_{\mathcal{E}}$ for each of the \mathcal{S}$$\mathcal{E} state-pairs within the microcanonical ensemble.

It is interesting that the derivation of $T_{\mathcal{SE}}$ in Eqs. 18-23 used only the standard temperature definition in Eqs. 12 and 18 and the choice of the zero-order basis for the densities of states $\rho_{\mathcal{E}}$ and $\rho_{\mathcal{SE}}$, used to formulate the sum in Eq. 19. In this respect the relation Eq. 23 is completely general, so it could also be used for other \mathcal{S}$$\mathcal{E} thermodynamic models which could potentially be much different from the simple oscillator model we use here.

III.3 Continuously varying time-dependent temperature

The temperature relations in the previous sections were derived using the standard expression Eq. 11 for the microcanonical ensemble, applied to the initial and final equilibrium states of the \mathcal{S}$$\mathcal{E} universe. It is useful to consider a time-dependent generalization of the microcanonical temperature that can be defined during thermalization. This uses time-dependent system probabilities from the system reduced density operator $\hat{\rho}_{\mathcal{S}}(t)$ in place of the microcanonical probabilities in Eq. 23, giving

[TABLE]

where $\rho_{\mathcal{S}}^{{n},{n}}$ is the probability of the system energy level $E_{n}$. Note that this time-dependent temperature agrees with the initial temperature $T_{\mathcal{E}}$ in Eq. 13 and with the final temperature $T_{\mathcal{SE}}$ in Eq. 23. $T_{\mathcal{SE}}(t)$ is the “master temperature” that describes the entire equilibration and thermalization process. Using Eq. 24 we will be able to follow the time-dependent changes in temperature as $\mathcal{S}$ and $\mathcal{E}$ begin in the initial state, exchange energy during thermalization, and eventually reach thermal equilibrium. This $T_{\mathcal{SE}}(t)$ is what we will be looking at as the “temperature” throughout the simulation.

IV Initial states for the simulations

The calculations start at $t=0$ with separable \mathcal{S}$$\mathcal{E} initial states

[TABLE]

where the initial system level is $|n_{0}\rangle$ and the initial environment state $|\epsilon_{0}\rangle$ has Gaussian distributed basis state probabilities

[TABLE]

with $\sigma_{\mathcal{E}}=0.5$ (the results are similar for other $0.1\leq\sigma_{\mathcal{E}}\leq 1$ that we have tested). In Eq. 26 the environment state is centered at an energy

[TABLE]

which varies with $n_{0}$, so that we are able to generate states that have the same nominal \mathcal{S}$$\mathcal{E} central energy $E_{0}=E_{\epsilon_{0}}+E_{n_{0}}$ but different system levels $n_{0}$. This will be useful for examining temperature equilibration, where the final state in principle will depend on the total energy but not on $n_{0}$. An example of the total probability per unit energy for an $n_{0}=4$ initial state $|\Psi_{n_{0}}\rangle$ at energy $E_{0}=5$ is shown in the top of Fig. 5. Each histogram bar in the figure shows the sum of \mathcal{S}$$\mathcal{E} basis states probabilities within the surrounding zero-order energy unit; the actual state is naturally much more complex in the zero-order basis. Note the logarithmic scale in the figure; the state is pretty sharply peaked around its nominal central energy. A slight asymmetry can be observed about the central energy $E_{0}=5$. This is because there are more basis states per unit energy above $E_{0}$ than below due to the increasing environment density of states. The asymmetry makes the average energy of the state slightly larger than the nominal energy $E_{0}$ in a way that depends on the environment density, which in turn depends on the environment energy $E_{\epsilon_{0}}$ and the system level $n_{0}$. This gives a slightly different initial state energy for each $n_{0}$, but the energies are close to the same.

We next consider the time evolution of this state, first with a random matrix coupling which we will find leads to pathological behavior, then with a more refined coupling that will be found to give physically satisfactory results.

V random matrix coupling and runaway thermalization dynamics

In this section we begin developing the quantum dynamics with the coupling Hamiltonian $\hat{H}_{\mathcal{SE}}$ of Eq. 1. We begin with a standard type of coupling, the random matrix coupling, used to model systems with classically chaotic dynamics Deutsch (2018), and often invoked in accounting for the existence of thermalization in quantum thermodynamics Esposito and Gaspard (2003); Borowski et al. (2003); Deutsch (2018). We used this in earlier simulations Barnes and Kellman (2013); Barnes et al. (2018); Lotshaw and Kellman (2019) with good results. However, we find here that with the introduction of a variable temperature, the random coupling introduces pathological behavior of runaway spreading of the wave packet. Furthermore, the random coupling is a serious limitation in itself – many important real systems do not have a random coupling. Thus, to understand thermalization for more realistic systems, we will want to explore more discriminating coupling forms.

The construction of $\hat{H}_{\mathcal{SE}}$ in Eq. 1 as a random matrix coupling begins with a matrix $\hat{R}$ filled with off-diagonal elements

[TABLE]

The $R_{n\epsilon,n^{\prime}\epsilon^{\prime}}$ are random complex numbers $R_{n\epsilon,n^{\prime}\epsilon^{\prime}}=X_{n\epsilon,n^{\prime}\epsilon^{\prime}}+iY_{n\epsilon,n^{\prime}\epsilon^{\prime}}$ as in Ref. Borowski et al. (2003). This is more generic than our previous work in Refs. Barnes and Kellman (2013); Barnes et al. (2018); Lotshaw and Kellman (2019), where we used real $R_{n\epsilon,n^{\prime}\epsilon^{\prime}}$ to minimize numerical effort. We generate the real and imaginary parts $X_{n\epsilon,n^{\prime}\epsilon^{\prime}}$ and $Y_{n\epsilon,n^{\prime}\epsilon^{\prime}}$ each as random numbers from a Gaussian distribution with standard deviation $\sigma=1$ with probabilities

[TABLE]

and similarly for the imaginary parts $Y_{n\epsilon,n^{\prime}\epsilon^{\prime}}$. We set the diagonal elements to zero to preserve the oscillator energies in the zero-order basis, as was done previously in Ref. Lotshaw and Kellman (2019). The interaction Hamiltonian is then constructed by multiplying $\hat{R}$ by a parameter $k$ that sets the overall coupling strength, $\hat{H}_{SE}=k\hat{R}$. This multiplication scales the random numbers so that their standard deviation becomes $\sigma=k$, consistent with the description in our earlier work Barnes and Kellman (2013); Barnes et al. (2018); Lotshaw and Kellman (2019) (e.g. in Eq. 10 of Ref. Barnes and Kellman (2013)). We chose $k$ to be the size of the average level spacing of the system-environment universe at our initial state energy $E_{0}=5$, since we have found that smaller $k$ do not give proper thermalization.

Fig. 5 shows time evolution with this coupling. With this coupling the initial Gaussian state associated with the top panel evolves in time to the state of the bottom panel. The time evolution evidently leads to runaway spreading of the wavepacket with probability in high energy states that does not appear to be converging to zero. This is not how a physically reasonable state should behave.

It is important to understand why this coupling causes runaway behavior here, because it was not observed, at least so prominently, in our earlier simulations with a fixed temperature bath. The coupling causes some spreading of the wavepacket to basis states of all energies, with the amount of probability per basis state decreasing rapidly as the states get farther off resonance from the initial state energy $E_{0}=5$. This might seem to entail decreasing probabilities at the top edge of the basis. However, the number of $\mathcal{E}$ basis states per unit energy increases very rapidly with increasing energy in the variable temperature bath, as shown in Fig. 1, so that many more basis states contribute to the total probability in each successive energy unit. Taken together, the total probability per unit energy doesn’t converge to zero as it should, as clearly seen in Fig. 5. This runaway coupling is a problem that needs to be addressed next.

VI Selective coupling “tames” thermalization dynamics

We will see that by defining a suitably much more selective coupling, physical results are obtained with both thermalization and contained spreading of the time-dependent quantum \mathcal{S}$$\mathcal{E} state. The basic idea is to “tame” the coupling to limit the range of transitions, especially to high energy states.

As before with the random matrix coupling, we begin with a coupling constant $k$ and a random matrix $\hat{R}$ as in Eq. 28. To construct $\hat{H}_{\mathcal{SE}}$, we take each individual matrix element of $k\hat{R}$ and multiply it by an exponential “taming” factor that depends on the quantum number differences between the coupled states:

[TABLE]

where $|\Delta n|=|n-n^{\prime}|$ is the quantum number difference between the coupled system states and $\sum_{osc}|\Delta n_{osc}|$ is the total quantum number difference for the individual oscillators in the coupled environment states. The parameters $\gamma_{\mathcal{S}}$ and $\gamma_{\mathcal{E}}$ suppress the coupling between \mathcal{S}$$\mathcal{E} states depending on how much they vary in quantum number, for example the coupling that moves one quantum between the system and bath is stronger than the coupling that moves two quanta. This limits the strength of transitions to high energy states, since they typically differ significantly in their quantum number distributions, thereby addressing the runaway problem.

A coupling scheme similar to Eq. 30 has been put forward by Gruebele Bigwood and Gruebele (1995); Gruebele (2003) in the context of intramolecular vibrational energy transfer, where he has argued that the exponential quantum-number dependence of the coupling is an approximate generic feature in molecular vibrational systems. Deutsch Deutsch (2018) has also said that a similar exponentially-tamed random matrix coupling can be obtained through a second-order perturbation theory analysis and that the exponential taming is needed to prevent runaway behavior in large quantum thermodynamic systems.

The tamed coupling has three parameters $k,\gamma_{\mathcal{S}},$ and $\gamma_{\mathcal{E}}$ that we choose somewhat arbitrarily for our model, with an aim towards obtaining physical thermalization behavior. The $k$ sets the “baseline” coupling strength; if $k$ is too small then thermalization will be impossible. The $\gamma_{\mathcal{E}}$ restricts the $\mathcal{E}$ transitions to address the runaway problem; it must be large enough to restrict the spreading with large energy differences, as needed for convergence, but also small enough to allow transfer between nearby $\mathcal{E}$ levels, as needed for thermalization. The $\gamma_{\mathcal{S}}$ controls how easily the system can transition between its levels; it must be small enough that all of the system levels can be accessed during the dynamics.

In our simulations we choose a coupling constant $k=0.15$. This is much larger than the $k$ we used with the random matrix coupling, to balance the exponential taming factors. We choose a relatively small system taming factor $\gamma_{\mathcal{S}}=0.125$ and a large environment factor $\gamma_{\mathcal{E}}=1$. This parameter choice gives good system thermalization behavior while limiting the environment transitions strongly enough to get good convergence within our basis. The effectiveness of this coupling and parameter choice is demonstrated by the time-evolved state in Fig. 6. The state corresponding to this figure began as an initial Gaussian state as seen in the top of Fig. 5, then it was evolved in time to equilibrium under the full Hamiltonian Eq. 1 containing the tamed coupling interaction $\hat{H}_{\mathcal{SE}}$ from Eq. 30. As seen in the histogram boxes in Fig. 6, the total probability per unit energy is converging to zero at the top edge of the basis. This shows that the tamed coupling has fixed the runaway problem of the random matrix coupling that was seen in the bottom of Fig. 5. Using the tamed coupling we found good convergence with a maximum \mathcal{S}$$\mathcal{E} energy $E^{max}=13$ for the simulations in this paper.

VII Results: Equilibration and Thermalization in the Simulations

Now we examine key aspects of the system dynamics during the approach to equilibrium: behavior of the time-dependent temperature; and the question of equilibrated Boltzmann distribution with thermalization. Is there thermodynamic-like behavior? But do we also see anomalous small-size temperature effects suggested by Fig. 3?

VII.1 Variable temperature and small-size effects

First we consider the computed time evolution of a set of initial states, constructed as described in Section IV with different initial system levels $n_{0}$ but the same nominal energies $E_{0}=6$. The total energies for the various $n_{0}$ are somewhat larger, as discussed in Section IV, with $6.116\leq\langle\hat{H}\rangle\leq 6.156$, where $\hat{H}$ is the total Hamiltonian Eq. 1. Taking $E=\langle\hat{H}\rangle$ in Eq. 23 we get for these states a narrow range of equilibrium microcanonical temperatures $1.912\leq T_{\mathcal{SE}}\leq 1.922$. Roughly speaking, we can think of all the states as sharing the common energy $E\approx 6.14$, hopefully corresponding in the simulations to a common final equilibrium temperature $T_{\mathcal{SE}}\approx 1.92$, where $1/T_{\mathcal{SE}}$ is the weighted average over all the initial state $1/T_{\mathcal{E}}$ at the common energy $E$, as in Eq. 23. We therefore test in the simulations whether the time-dependent temperature $T_{\mathcal{SE}}(t)$ of Eq. 24 equilibrates to the common temperature $T_{\mathcal{SE}}\approx 1.92$.

Fig. 7 shows the time-dependent behavior of the temperatures $T_{\mathcal{SE}}(t)$ for each of the initial states $n_{0}$. For each $n_{0}$, the temperature begins in its respective value for an isolated system and environment, $T_{\mathcal{SE}}(t=0)=T_{\mathcal{E}}$ (from Eqs. 24 and 13). Time evolution takes the temperatures to equilibrium, where they do in fact fluctuate around the common approximate value $T_{\mathcal{SE}}\approx 1.92$. Thus, we are getting the common microcanonical $T_{\mathcal{SE}}$ value corresponding to energy $E\approx 6.14$, as hoped for. This result validates the path of development in Section III regarding a variable temperature. Observed small temperature fluctuations at equilibrium are due to the time-dependent fluctuations in the system density operator $\hat{\rho}_{\mathcal{S}}(t)$, whose behavior will be discussed shortly in Section VII.2.

It is a noteworthy prediction based on the considerations of Section III that the finite bath equilibrium temperatures in Fig. 7 should be considerably higher than would be expected using the infinite bath $T$ from Eq. 14 based on the average number of quanta per degenerate oscillator $\langle n_{osc}\rangle=$ $\langle E_{osc}\rangle$. To test this, we calculated $\langle E_{osc}\rangle=\langle E_{\mathcal{E}}\rangle/\eta$ as the time-averaged equilibrium value for times $30<t\leq 60$ averaged over all of the simulations shown in Fig. 7, giving $\langle E_{osc}\rangle=1.117\pm 0.004$. The infinite bath limit temperature Eq. 14 from this $\langle E_{osc}\rangle$ is $T=1.564\pm 0.004$, much smaller than our temperature $T_{\mathcal{SE}}=1.92$. This is because the finite bath temperatures $T_{\mathcal{E}}$ in Eq. 13 (which go into the calculation of the $T_{\mathcal{SE}}$ via Eq. 23) increase more rapidly with energy than the infinite bath $T,$ as was seen in Fig. 3. Thus, the anomalous temperature scaling of the small environment is demonstrably evident from this analysis of Fig. 7. We will have more to say about the anomalous temperature in the next subsection.

VII.2 Approach to thermal equilibrium and anomalous size effects

Next, we consider the behavior of the system in the approach to thermal equilibrium. Fig. 8 shows an example of the time-dependent system probabilities $\rho_{\mathcal{S}}^{n,n}$ from the reduced density operator for an initial $\mathcal{S}$ level $n_{0}=0$ (the dynamics are similar for the other $n_{0}$). As the state begins to evolve in time, much of the initial state probability is quickly lost to the other levels, followed by a much slower decay to the equilibrium Boltzmann distribution marked by the dotted lines. The behavior can be fit by an empirical power law

[TABLE]

where $\tau$ and $\delta$ are fit parameters and $\exp(-E_{n_{0}}/T_{\mathcal{SE}})/Z$ is the equilibrium Boltzmann probability at the temperature $T_{\mathcal{SE}}$, as will be discussed further shortly. Power law decays have been discussed by Gruebele Gruebele (1998, 2003) as a generic feature in molecular vibrational systems that can be described by couplings similar to our Eq. 30. The decay describes the nearly exponential drop of the initial state $n_{0}$ probability at short times and the longer decay to equilibrium. The other levels $n$ reach equilibrium at different timescales depending on how far they are from the initial level $n_{0}=0$, for example, $n=1$ reaches its equilibrium probability relatively quickly whereas it takes much longer for the $n=4$ level. This stands in contrast to the dynamics under the simple random matrix coupling, where each system level evolves at approximately the same rate Barnes and Kellman (2013), without any sense of “proximity” between nearby energy levels that facilitates their energy transfer. Beyond simply being essential to converge the calculations, as discussed in Section VI, it seems to us that the tamed coupling is also giving a much more realistic dynamics .

At long times, the system level probabilities fluctuate about a Boltzmann-appearing distribution $\rho_{\mathcal{S}}^{n,n}$ $\sim\exp(-E_{n}/T_{\mathcal{SE}})$ at the temperature $T_{\mathcal{SE}}$, shown as a black dotted line for each $E_{n}$. The agreement with the Boltzmann distribution at $T_{\mathcal{SE}}$ is examined in Fig. 9 across a range of initial state energies $E=\langle\hat{H}\rangle$ and corresponding temperatures listed in Table 2. The time-averaged system probabilities from the simulations are in very good agreement with the analytical Boltzmann distributions at temperatures $T_{\mathcal{SE}}$ from Eq. 23. For comparison, in Fig. 9 we also show the Boltzmann distributions for the infinite bath temperatures $T$ calculated for the states, based on the average energy per bath oscillator observed in the simulations, see Table 2 and the discussion in the last paragraph of Section VII.1. The resulting temperatures are systematically lower than the $T_{\mathcal{SE}}$ values, and the corresponding Boltzmann distributions do a poor job of describing the system probabilities. Thus, the observed thermalization to $T_{\mathcal{SE}}$ strongly reinforces that this is the correct thermodynamic temperature to describe the total system \mathcal{S}$$\mathcal{E} .

At this point it is appropriate to remark on the question of “eigenstate thermalization” in our simulations. The eigenstate thermalization hypothesis (ETH), that eigenstates of a suitable system-environment Hamiltonian reflect thermal properties Deutsch (2018, 1991); D’Alessio et al. (2016), is widely regarded as an explanation for thermalization phenomena. ETH is often justified through an appeal to chaotic dynamics of the kind that classically corresponds to a random matrix Hamiltonian. Chaotic dynamics become less certain the more that there is a “taming” of the coupling, as used in this paper to get convergence of the dynamics, and ETH thereby becomes less certain as well. Nonetheless, all of our initial states thermalize to their expected temperatures, and this is consistent with ETH. In future work, we plan to explore the breakdown of ETH as reduced coupling strength makes questionable chaotic dynamics, ETH behavior, and thermalization itself.

Another point worth remark is alternatives to the random matrix-based couplings used in this paper. Simple couplings based on linear combinations of raising and lowering operators are used in many quantum thermodynamic investigations D’Alessio et al. (2016). Accordingly, we have run calculations where we adopt a linear $k\hat{x}_{i}\hat{x}_{j}$ coupling. We find that this gives controlled spreading with semi-quantitative thermalization. However, in comparison the thermalization is significantly better with the random matrix tamed coupling calculations reported above. The likely reason the random matrix works better for our setup is that our five-oscillator bath has approximate frequency resonances. This is typical of many physical systems, e.g. a molecule embedded in a bath, which will almost inevitably have such “anharmonic resonances.” A random coupling will better capture the effects of these resonances. On the other hand, there are systems, e.g. of coupled bosons, where the $\hat{x}_{i}\hat{x}_{j}$ type coupling is more appropriate. Based on our calculations, we believe that variable temperature baths can be devised appropriate to a variety of physical situations in “tailor-made” fashion.

VIII summary and prospects

This paper has considered a quantum description of energy flow from a system into a very small variable temperature bath. We defined a system, consisting of a finite number of levels, and an environment, consisting of levels of a finite collection of harmonic oscillators (which constitutes the bath). A set of identical oscillators was first considered, paralleling the Einstein heat capacity model. To get something more like a continuous state distribution, we then took a collection of non-identical oscillators. This gives a distribution of levels that closely tracks that of the bath of identical oscillators, but also has the desired feature of breaking the degeneracy, giving a quasi-continuous level distribution. The level pattern of this bath has a density of states that gives temperature-like behavior, using the standard statistical thermodynamic microcanonical relation between temperature, energy, and density of states. This defines the “temperature” $T_{\mathcal{E}}$ for the finite bath. This temperature differs significantly from that of the infinite oscillator bath, as seen in simulations with a bath with only $\eta=5$ oscillators. We compared the energy-temperature relations for a single oscillator within the infinite bath (the well-known result of Einstein from his famous heat capacity paper) to the corresponding relation for a finite bath. There are systematic differences, which are pronounced for $\eta=5$, and asymptotically approach the infinite bath at large $\eta$. The small bath has higher temperature for a given amount of energy per oscillator. Very unlike the infinite bath, it also terminates at a temperature $T_{\mathcal{E}}>0$, as seen in Fig. 3.

Having devised the finite bath with temperature $T_{\mathcal{E}}$, we considered the process of heat flow from the system into this bath. Simulations were performed of the process of heat flow to the finite bath in quantum time evolution. First we used a random-matrix coupling of the kind that has been employed in many contexts, including successful quantum thermodynamic simulations Borowski et al. (2003); Barnes and Kellman (2013); Barnes et al. (2018); Lotshaw and Kellman (2019). This however led to “runaway spreading” of the quantum \mathcal{S}$$\mathcal{E} wave function. This is closely connected with the variable temperature of the bath – a feature not present in earlier thermodynamic simulations. The problem is that the density of states increases rapidly with increasing temperature, and the non-discriminate random coupling overpowers the quantum time evolution. To solve this, we switched to a more selective coupling similar to the kind that has long been used Bigwood and Gruebele (1995); Gruebele (2003) in molecular simulations. This selective coupling “tames” the spreading of the wave function, so that runaway behavior is avoided. The tamed coupling appears to be a realistic new feature needed to solve a real problem in the simulations.

Next came computational examination of the temperature $T_{\mathcal{SE}}$ defined for the microcanonical ensemble of the \mathcal{S}$$\mathcal{E} total system “universe,” including the time-dependent temperature $T_{\mathcal{SE}}(t)$ that varies continuously between the initial bath temperature $T_{\mathcal{E}}$ and the final \mathcal{S}$$\mathcal{E} temperature $T_{\mathcal{SE}}$. In simulations with the $\eta=5$ oscillator bath, starting with different initial system states but the same total system-environment energy, we tracked the temperature from its various initial values (because the bath has different energies depending on the system state) to its final value at equilibrium. All the simulations went to essentially the same final temperature $T_{\mathcal{SE}}$, as desired. The simulations with the bath of $\eta=5$ oscillators with selective coupling show equilibration to a Boltzmann-type distribution at the temperature $T_{\mathcal{SE}}$ implied by the initial energy of the total system. As noted above, this temperature is markedly different from that of an infinite bath with the equivalent energy per bath oscillator. In short, there are marked effects of the small finite bath on thermal behavior with variable temperature in the quantum simulations.

It is interesting to consider real situations in which to explore these finite size quantum thermodynamic effects. Experiments on very small Bose-Einstein condensates, containing as few as six atoms Kaufman et al. (2016), may point the way to size-dependent variable temperature behavior similar to the oscillator model we have studied here. Several investigators have proposed small molecules as laboratories for fundamental exploration of quantum thermodynamics and statistical mechanics. Leitner Leitner (2015, 2018) has reviewed a method of using the eigenstate thermalization hypothesis to understand ergodicity and localization of energy within time-dependent molecular systems. Pérez and Arce Pérez and Arce (2018) performed simulations of dynamics on a potential energy surface of the molecule OCS, which has a long history as an exemplar of problems of classically chaotic molecular dynamics. They treat one of the vibrational modes of OCS as a “system,” and the other two modes as an “environment,” akin to what we do here, but with a two-mode bath that is much smaller even than what has been considered here. They find a kind of thermalization of the system when it is excited with sufficient energy to have chaotic classical dynamics. However, they did not engage in the kind of analytic treatment of temperature of the present paper. If we go to a four-atom molecule, for example the important species C2H2 (acetylene) or H3O+ (hydronium ion), we could take as system one of the modes, e.g. a C-H stretch, leaving 5 vibrational modes as the bath, just as we do here. This ignores rotational degrees of freedom; one could do experiments with angular momentum $J$ = 0; or alternately, allow $J$ excitations, which would become increasingly important at higher $J$, where rotation-vibration coupling would become important, giving the rotational degrees of freedom as a second bath or environment $\mathcal{E}$′. It is worth noting that molecular systems interacting with small baths are of interest in other contexts as well, e.g. in calculations of entanglement dynamics and spectroscopic signals Cheng and Cina (2014); Kovac and Cina (2017).

As an alternative to the molecular dynamics simulations of Ref. Pérez and Arce (2018), one could also use “effective Hamiltonians” of the kind that have had vast use in molecular spectroscopy Kellman (1995); Tyng and Kellman (2007). It is notable that these Hamiltonians usually employ one or more “polyad numbers” that constitute approximate constants of motion, valid on a limited time scale. This makes these attractive systems in which to explore the effects of approximate constants as barriers to thermalization, a topic of considerable interest Deutsch (2018) in contemporary theory of quantum thermodynamics. The effective molecular polyad Hamiltonian can then be enhanced with polyad-breaking perturbations Chakraborty and Kellman (2008); Barnes and Kellman (2010, 2011) that correspond to real molecular dynamical effects. These hierarchical dynamical systems could be ideal laboratories for investigation of thermodynamic processes on multiple time scales.

As a final comment, taking a wider perspective on the work here, it may be worthwhile to consider that there are (at least) three dimensions of “post-classical” effects in quantum thermodynamics. The first of course is quantization of energy levels, introduced in the very beginnings of quantum physics by Planck in his black-body theory and by Einstein in his famous heat capacity paper. A second is finite size, as exemplified in this paper by the very small size (five oscillators) of the variable temperature bath. A third involves quantum time evolution. This might come with more complicated setups of finite size and time evolution than explored here. One might consider a system linking two baths of different sizes; or a system linking two finite baths where the coupling of the system to each bath is different. These would require far larger simulations than performed here. We can readily imagine experimental realizations of these situations, e.g. with supramolecular arrangements of two or more molecules weakly linked by a third.

IX Acknowledgements

P. L. thanks Rob Yelle and Craig Rasmussen for technical assistance with computations. This work benefited from access to the University of Oregon high performance computer Talapas.