Quantum Joule Expansion of One-Dimensional Systems

TL;DR

This paper studies the quantum Joule expansion in one-dimensional nonintegrable systems, revealing thermalization, negative effective temperatures, and the potential to create negative temperature states dynamically, with implications for quantum thermodynamics.

Contribution

It demonstrates how quantum Joule expansion leads to thermalization and negative temperatures in 1D systems, and proposes a method to generate negative temperature states dynamically.

Findings

Diagonal and canonical ensemble weights match at high initial temperatures.

Subsystem entropies equilibrate to thermal entropy at long times.

Negative effective temperatures occur in finite systems and vanish in the thermodynamic limit for bosons.

Abstract

We investigate the Joule expansion of nonintegrable quantum systems that contain bosons or spinless fermions in one-dimensional lattices. A barrier initially confines the particles to be in half of the system in a thermal state described by the canonical ensemble and is removed at time $t = 0$. We investigate the properties of the time-evolved density matrix, the diagonal ensemble density matrix and the corresponding canonical ensemble density matrix with an effective temperature determined by the total energy conservation using exact diagonalization. The weights for the diagonal ensemble and the canonical ensemble match well for high initial temperatures that correspond to negative effective final temperatures after the expansion. At long times after the barrier is removed, the time-evolved R\'enyi entropy of subsystems bigger than half can equilibrate to the thermal entropy with exponentially small fluctuations. The time-evolved reduced density matrix at long times can be approximated by a thermal density matrix for small subsystems. Few-body observables, like the momentum distribution function, can be approximated by a thermal expectation of the canonical ensemble with strongly suppressed fluctuations. The negative effective temperatures for finite systems go to nonnegative temperatures in the thermodynamic limit for bosons, but is a true thermodynamic effect for fermions, which is confirmed by finite temperature density matrix renormalization group calculations. We propose the Joule expansion as a way to dynamically create negative temperature states for fermion systems with repulsive interactions.

Figures

Equations

$${\hat{H}}^b=-J\sum _{l=1}^{L-1}(a_l^{†}{a}_{l+1}+H.c.)+\frac{U}{2}\sum _{l=1}^L{n}_l^b(n_l^b-1),$$

$${\hat{H}}^f$$

$${\hat{\rho }}^0=\sum _i{\rho }_i^0|i⟩⟨i|$$

$${\hat{\rho }}^0$$

$${\hat{\rho }}^t=\sum _{m,n}{\rho }_{mn}^0\exp [-it(E_m-E_n)]|m⟩⟨n|.$$

$$\hat{\bar{\rho }}={\hat{\rho }}^d$$

$$⟨\hat{O}{⟩}_t$$

$$⟨{\hat{H}}_F⟩$$

$$\frac{dE(\beta )}{d\beta }=\sum _{m,n}-\frac{e^{-\beta (E_n+E_m)}(E_n-E_m{)}^2}{2Z^2}<0,$$

$${\hat{\rho }}^c=\sum _m\frac{e^{-{\beta }_{eff}{E}_m}}{Z({\beta }_{eff})}|m⟩⟨m|=\sum _m{W}_m^c|m⟩⟨m|.$$

$$⟨\hat{O}⟩={\text{Tr}}_A\left(\hat{O}{\hat{\rho }}_A\right).$$

$$S_{pA}=\frac{1}{1-p}\ln {\text{Tr}}_A{\left({\hat{\rho }}_A\right)}^p.$$

$$S_{2A}^t=$$

$$\Delta E$$

$$S_{2A}^t\approx$$

$$-\ln (\sum _{m,m^{\prime },i_B,i_B^{\prime }}{\rho }_{mm}^0{\rho }_{m^{\prime }{m}^{\prime }}^0⟨m|i_B⟩⟨i_B^{\prime }|m^{\prime }⟩⟨m^{\prime }|i_B^{\prime }⟩⟨i_B|m⟩$$

$$+\sum _{m\ne m^{\prime },i_B,i_B^{\prime }}{\rho }_{m{m}^{\prime }}^0{\rho }_{m^{\prime }m}^0⟨m^{\prime }|i_B⟩⟨i_B^{\prime }|m^{\prime }⟩⟨m|i_B^{\prime }⟩⟨i_B|m⟩)$$

$$=-\ln (\sum _{m,m^{\prime }}{\rho }_{mm}^0{\rho }_{m^{\prime }{m}^{\prime }}^0{T}_{AB}^{m{m}^{\prime }}$$

$$+\sum _{m\ne m^{\prime }}{\rho }_{m{m}^{\prime }}^0{\rho }_{m^{\prime }m}^0{T}_{BA}^{m^{\prime }m})$$

$$T_{BA}^o\sum _{m\ne m^{\prime }}{\rho }_{m{m}^{\prime }}^0{\rho }_{m^{\prime }m}^0=T_{BA}^o\left(e^{-S_2^{{\beta }_0}}-e^{-S_2^d}\right)$$

$$S_{2A}^t\approx S_{2A}^d-T_{BA}^o{e}^{S_{2A}^d-S_2^{{\beta }_0}}$$

$$Var(S_{2A}^t)$$

$${\hat{\rho }}^c(AB)=$$

$${\hat{\rho }}^c(A)$$

$$n_k^f=\frac{1}{L}\sum _{p,q}{e}^{-ik(p-q)}⟨c_p^{†}{c}_q⟩,$$

$$n_k^b=\frac{1}{L}\sum _{p,q}{e}^{-ik(p-q)}⟨a_p^{†}{a}_q⟩.$$

$${△}^{2}O$$

$${△}^{2}O<M\left(e^{-S_2^{{\beta }_0}}-e^{-S_2^d}\right)$$

Full text

I Introduction

With the remarkable advances in efficient computing algorithms and cold atom experiments, nonequilibrium dynamics has been extensively studied both theoretically and experimentally in recent years. Quantum thermalization is one of the most important topics in this area. In 1929, the classical theory of statistical mechanics was reformulated quantum-mechanically by von Neumann von Neumann (2010), which opened the door to the study of quantum thermalization through the unitary dynamics of quantum systems. Later studies Jensen and Shankar (1985); Deutsch (1991); Goldstein et al. (2010) show that quantum thermalization is closely related to random matrix theory Mehta (2004); Guhr et al. (1998) and the eigenstate thermalization hypothesis (ETH) Srednicki (1999); Rigol and Srednicki (2012) was proposed to understand it in finite quantum systems. For nonintegrable quantum systems with Wigner-Dyson level spacing distribution Santos and Rigol (2010a), ETH asserts that every eigenstate is thermal and, as a result, few-body observables thermalize Kim et al. (2014); Yoshizawa et al. (2018). ETH has been verified for a variety of nonintegrable quantum systems Rigol et al. (2008); Rigol (2009); Santos and Rigol (2010b); Biroli et al. (2010); Genway et al. (2012); Khatami et al. (2012); Neuenhahn and Marquardt (2012); Khatami et al. (2013); Ikeda et al. (2013); Kim et al. (2014); Beugeling et al. (2014); Sorg et al. (2014); Beugeling et al. (2015); Mondaini et al. (2016); Mondaini and Rigol (2017); Yoshizawa et al. (2018); Hikida et al. (2018), with exceptions when little entanglement in the eigenbasis results in equilibration without thermalizationGogolin et al. (2011). A modified version of ETH for subdiffusive thermalization Luitz and Bar Lev (2016) and strong forms of ETH for reduced density matrix Garrison and Grover (2018); Dymarsky et al. (2018) have been formulated recently.

The quantum quench protocol is often used in studies of nonequilibrium dynamics, where an initial state of the system is prepared at time $t=0$, some parameters are then suddenly or time-dependently changed and the state is let to evolve under the new Hamiltonian. This operation involves a global quench Calabrese and Cardy (2006, 2007a); Sotiriadis et al. (2009); Rigol (2009); Rigol and Fitzpatrick (2011); He and Rigol (2012, 2013); Sorg et al. (2014); Mistakidis et al. (2014, 2015); Cardy (2016); Neuhaus-Steinmetz et al. (2017); Mistakidis and Schmelcher (2017); Plaßmann et al. (2018); Mistakidis et al. (2018), local quench Calabrese and Cardy (2007b); Eisler and Peschel (2007) or sudden expansion Rigol and Muramatsu (2004); Minguzzi and Gangardt (2005); Rigol and Muramatsu (2005a, b); Camalet (2008); Heidrich-Meisner et al. (2009); Kajala et al. (2011); Vidmar et al. (2013); Ronzheimer et al. (2013); Sorg et al. (2014); Xia et al. (2015); Campbell et al. (2015); Vidmar et al. (2015, 2017a, 2017b); Xu and Rigol (2017); Herbrych et al. (2017); Koutentakis et al. (2017); Noh et al. (2018); Scherg et al. (2018); Siegl et al. (2018). The initial states can be either pure states or any mixed states. The dynamics of sudden expansion has rich results including dynamical fermionization in the expansion of a Tonks-Girardeau gas Minguzzi and Gangardt (2005); Vidmar et al. (2017b), hard-core bosons Rigol and Muramatsu (2005a, b); Xu and Rigol (2017) or a Bose gas Campbell et al. (2015), quasi-condensation at finite momentum in hard-core bosons Rigol and Muramatsu (2004, 2005a, 2005b); Vidmar et al. (2013, 2015); Xu and Rigol (2017), ballistic and diffusive expansion for bosons and fermions in one and two dimensions Kajala et al. (2011); Ronzheimer et al. (2013); Vidmar et al. (2013), quantum distillation of singlons and doublons during the expansion Heidrich-Meisner et al. (2009); Xia et al. (2015); Herbrych et al. (2017); Scherg et al. (2018), expansions of atomic clouds and interwell tunneling dynamics of a Bose-Fermi mixture Siegl et al. (2018) and self-trapping of bosons initially confined in an harmonic trap Koutentakis et al. (2017). A time-dependent emergent local Hamiltonian can be constructed to analytically describe these nonequilibrium states as its eigenstates Vidmar et al. (2017a, b); Xu and Rigol (2017).

A celebrated experiment in the context of classical statistical mechanics is the Joule expansion. The free expansion of an ideal gas from an initial volume $V$ to a final volume $2V$ does not change the temperature of the gas and the increase in entropy is $nR\ln{2}$. For interacting gases, the temperature decreases for attractive interactions, such as for the van der Waals gas, and increases for repulsive interactions. But what happens for quantum systems? The Joule expansion of an isolated perfect quantum gas is discussed in Camalet (2008), where the time evolution of the particle number density displays periodicity with period proportional to the square of the size of the system, consistent with the recurrence time for free bosons and free fermions obtained in Kaminishi et al. (2015), and much smaller than the classical Poincaré recurrence time, which scales exponentially with the size of the system. For generic quantum systems, the recurrence time Bocchieri and Loinger (1957) is much longer, typically scaling as a double exponential of the system size Kaminishi et al. (2015); Mori et al. (2018).

Sudden expansion of an initially confined thermal state for hard-core bosons was studied in Ref. Xu and Rigol (2017), where the dynamical fermionization and the quasi-condensation during the expansion are observed. These studies focus on integrable systems where analytic tools can be used. But discussions of the question of thermalization at long times are still missing, which is a prerequisite to answer the Joule question in the quantum regime.

In this article, we discuss the Joule expansion for two interacting but nonintegrable one-dimensional quantum Hamiltonians. We use thermal initial states for a range of temperatures. We discuss the thermalization after a sudden expansion where the particles are kept confined within a length which is twice the initial one. We discuss two important problems which to the best of our knowledge have not been discussed in the context of sudden expansion in the ETH literature. The first one is to establish if the long-time evolution can be characterized by the diagonal ensemble or the canonical ensemble for a final temperature and to which extent these ensemble descriptions are compatible. The second one is the long-time behavior of the Rényi entropies for subsystems of various sizes. If the subsystem is large enough, can these entropies stabilize to thermal values? In addition, we will discuss the momentum distributions before and after the expansion as these can be measured experimentally by time-of-flight experiments with cold atoms trapped into one-dimensional optical lattices often called “tubes”.

The paper is arranged as follows. In Section II we introduce the boson and spinless fermion models that we study and discuss the time evolution following the removal of the barrier (expansion from size $L/2$ to size $L$), the diagonal and the canonical ensemble descriptions for observables in the equilibrated regime, and the Rényi entropy. In Section III, the weights of the eigenstates for the diagonal ensemble (DE) and the canonical ensemble (CE) descriptions are calculated and compared. We find that, above a certain initial temperature, the final effective temperature obtained from the canonical ensemble description is negative. Results for the von Neumann entropy and the second order Rényi entropy are presented in Section IV. In Section V we investigate the full long-time evolved solution, the diagonal ensemble description and the canonical ensemble with an effective final temperature at the level of the reduced density matrices and their eigenvalues and eigenstates in these three cases. In Section VI, we calculate the momentum distribution function, which at long times shows the expected population inversion for the negative effective temperature cases. In Section VII we present the answer to the Joule question, that is, the final effective temperature after expansion from $L/2$ to $L$ as function of the temperature of the thermal state at $t=0$. Finite-size analysis shows that effective negative temperatures survive in the thermodynamic limit for the case of spinless fermions, suggesting that Joule expansion can be used as a way to dynamically create negative temperature fermionic states. Finally, Section VIII contains a summary and conclusions.

II Model and Methods

We consider two models: the one-dimensional Bose-Hubbard model and one-dimensional spinless fermions with nearest-neighbor (NN) and next-nearest-neighbor (NNN) hopping and interaction. Open boundary conditions are used for both models. The Hamiltonian for the Bose Hubbard model is:

[TABLE]

and the Hamiltonian for spinless fermions is:

[TABLE]

where $L$ is the total number of sites, $a^{\dagger}_{l}$ ($a_{l}$) is the creation (annihilation) operator of bosons at site $l$, satisfying the commutation relations $[a_{l},a^{\dagger}_{l^{\prime}}]=\delta_{ll^{\prime}}$, and $c^{\dagger}_{l}$ ($c_{l}$) is the creation (annihilation) operator of spinless fermions at site $l$, satisfying the anitcommutation relations $\{c_{l},c^{\dagger}_{l^{\prime}}\}=\delta_{ll\prime}$. $n^{b}_{l}=a^{\dagger}_{l}a_{l}$ ($n^{f}_{l}=c^{\dagger}_{l}c_{l}$) is the particle number operators for bosons (fermions) at site $l$. Repulsive interactions $U,V_{1},V_{2}>0$ are considered for both systems. For nonintegrable systems, the energy levels are repulsive and the level spacing distribution is given by the Wigner-Dyson distribution, which has been found for both $\hat{H}^{b}$ Kollath et al. (2010) and $\hat{H}^{f}$ Santos and Rigol (2010a). The fermionic model $\hat{H}^{f}$ contains next-nearest-neighbor terms $J_{2}$ and $V_{2}$ that may seem prohibitively difficult to achieve experimentally. We note, however, that this model can also be realized with a zigzag ladder, with one leg containing the odd sites of Eq. (2) and the other leg containing the even sites. The $J_{1}$ and $V_{1}$ in Eq. (2) are terms connecting odd and even sites, and therefore are the inter-leg zig-zag terms in the ladder system. In the ladder system, the $J_{2}$ and $V_{2}$ terms then become the intra-leg nearest-neighbor terms. We set $\hbar=k_{B}=1$ in the following and fix the hopping energy $J_{1}=V_{1}=J=1$ in our calculations. The results we present are for spinless fermions with $J_{2}=V_{2}=1$ and bosons with $U=3$, unless otherwise specified.

In the Joule expansion, a volume of gas is initially prepared with a thermal state inside one half of a container via a partition, with the other half of the container empty. The partition is then suddenly removed and the gas is let to expand freely to occupy the whole volume. For particles on a lattice, a similar initial setup can be achieved by adding a high potential throughout the right half of the lattice, $\hat{H}_{W}=W\sum_{l=L/2+1}^{L}n_{l}$, where $W\gg J,U,J_{1},J_{2},V_{1},V_{2}$. In the numerical calculations presented here, we set $W=10^{6}$ to make sure that $\beta W\gg 1$ in all calculations, so particles won’t hop to the right half due to high temperatures. At time $t=0$, the high potential is suddenly removed ($W$ is quenched to zero) and particles can then move on a twice bigger lattice. For fermions (bosons), we consider the initial state to be a $N_{p}=L/4$-particle ($L=20$ unless otherwise specified) thermal state $\hat{\rho}^{f}_{0}$ ($\hat{\rho}^{b}_{0}$) at a finite temperature $T^{f}_{0}=1/\beta^{f}_{0}$ ($T^{b}_{0}=1/\beta^{b}_{0}$) in the left half of the lattice.

The general language describing the time evolution of a quantum system is the following. Let $\{|i\rangle\}$ be the eigenstates of the initial Hamiltonian before expansion, $\hat{H}_{0}|i\rangle=e_{i}|i\rangle$, where $\hat{H}_{0}=\hat{H}^{b(f)}+\hat{H}_{W}$ for bosons (fermions). The initial mixed thermal state $\hat{\rho}_{0}$ with inverse temperature $\beta_{0}$ in this eigenbasis is given by

[TABLE]

where $\rho^{0}_{i}=e^{-\beta_{0}e_{i}}/Z_{0}(\beta_{0})$, and $Z_{0}(\beta_{0})=\textrm{Tr}{(e^{-\beta_{0}\hat{H}_{0}})}$ is the partition function.

Let $\{|m\rangle\}$ be the eigenstates of the final Hamiltonian $\hat{H}_{F}$ that drives the time evolution after quenching (removing) the high wall, that is, $\hat{H}_{F}|m\rangle=E_{m}|m\rangle$, with $\hat{H}_{F}=\hat{H}^{b(f)}$ for bosons (fermions). The initial state, written in this eigenbasis is given by,

[TABLE]

and contains diagonal ($m=n$) and off-diagonal ($m\neq n$) terms. The time-dependent density operator at time $t$ during the expansion is

[TABLE]

If there is no degeneracy in the spectrum, which is generally true for nonintegrable systems Bohigas et al. (1984), the long-time average of the density operator above will give the diagonal ensemble (DE) Rigol et al. (2008) density matrix

[TABLE]

where $\rho^{0}_{mm}$ is the weight $W^{d}_{m}$ that the projector $\hat{P}_{m}=|m\rangle\langle m|$ has in the DE.

For an arbitrary observable $\hat{O}=O_{m^{\prime}n^{\prime}}|m^{\prime}\rangle\langle n^{\prime}|$, the time-dependent expectation value is,

[TABLE]

where the first term $\langle\hat{O}\rangle_{d}=\sum_{m}O_{mm}\rho^{0}_{mm}$ is the expectation value in the diagonal ensemble, or the long-time average value. The second term, for nonintegrable systems, is expected to be very small at long times for two reasons: firstly, according to ETH, the off-diagonal matrix elements for a few-body observable $\hat{O}$ is exponentially small in the size of the system Deutsch (1991); Srednicki (1999), secondly the long time temporal dephasing Rigol et al. (2008) guarantees the canceling of oscillations at long enough times. It seems problematic at first that one may need to wait an exponentially long time for the cancellation due to dephasing to occur because of the exponentially small level spacings. It was however pointed out D’Alessio et al. (2016) that since the off-diagonal matrix elements of $\hat{O}$ are exponentially small, the phase coherence between only a finite fraction of the eigenstates with a significant contribution to the expectation value needs to be destroyed. If the observable $\hat{O}$ thermalizes, the expectation in the diagonal ensemble must agree with the expectation value obtained in statistical mechanics. We use the canonical ensemble (CE) to describe this thermalized state. We find the effective inverse temperature $\beta_{eff}$ of this CE by matching the total energy of the system, which is a time-independent conserved quantity,

[TABLE]

where $Z(\beta_{eff})=\textrm{Tr}{(e^{-\beta_{eff}\hat{H}_{F}})}$ is the partition function for the CE at temperature $T_{eff}=1/\beta_{eff}$. The function $E(\beta)=\textrm{Tr}{(\hat{H}e^{-\beta\hat{H}})}/Z(\beta)$ is monotonically decreasing with $\beta$ because

[TABLE]

and therefore there is only one solution for the effective temperature. For the CE we then have

[TABLE]

Like in a classical gas, the average energy of the initial Hamiltonian $H_{0}$ (particles confined to $L/2$) with repulsive interactions at a certain temperature is higher than that of the final Hamiltonian $H_{F}$ (system size $L$) with the same temperature. So a higher effective temperature is needed for Eq. (8) to be satisfied. Conversely, for attractive interactions (such as in the classical van der Waals gas), the system is expected to cool down upon expansion. For repulsive interactions, such as the systems studied here, an interesting possibility can occur. If the initial average energy is higher than the infinite-temperature energy for the final Hamiltonian, a negative temperature is needed to compensate for the difference. It is convenient to describe this in terms of inverse temperature: the system is prepared in some positive inverse temperature that decreases upon expansion for repulsive interactions. There will be some initial positive inverse temperature for which it decreases to zero (infinite temperature) after expansion. Any positive initial inverse temperature that is smaller than that will then lead to even smaller final inverse temperatures, leading to negative values.

We have the initial thermal state $\hat{\rho}^{0}$, the time-evolved state $\hat{\rho}^{t}$, the DE density matrix $\hat{\rho}^{d}$ and the corresponding CE density matrix $\hat{\rho}^{c}$. To better compare these states, we investigate their reduced density matrices. Note that for an arbitrary density operator $\hat{\rho}$, if the observable $\hat{O}$ only resides on a subsystem $A$ with a size $l_{A}$, whose complement is denoted as $B$ with a size $l_{B}$, the expectation value only depends on the reduced density operator of the subsystem $\hat{\rho}_{A}=\textrm{Tr}_{B}{\hat{\rho}}$ by

[TABLE]

If the reduced density matrix is similar to the CE reduced density matrix, every observable that resides in subsystem $A$ has the same expectation value as in the canonical ensemble.

Another characterization of the reduced density matrix is the Rényi entropy. The $p$-th order bipartite Rényi entropy for $\hat{\rho}_{A}$ is defined as

[TABLE]

The von Neumann entanglement entropy $S_{1A}=-\textrm{Tr}{\left(\hat{\rho}_{A}\ln{\hat{\rho}_{A}}\right)}$ is obtained by taking the limit $p\rightarrow 1^{+}$. For pure states, the Rényi entropy probes the entanglement between subsystems $A$ and $B$, and $S_{pA}=S_{pB}$. The conformal field theory (CFT) prediction of Rényi entropy for the ground state of critical one-dimensional systems is discussed in Holzhey et al. (1994); Vidal et al. (2003); Jin and Korepin (2004); Calabrese and Cardy (2004, 2009). For excited states of chaotic many-body Hamiltonians that have finite energy density in the sense that $\lim_{V\rightarrow\infty}(E_{n}-E_{0})/V\neq 0$ Garrison and Grover (2018), Ref. Lu and Grover (2019) argues that the von Neumann entanglement entropy is proportional to the subsystem size for $l_{A}/L<1/2$, while $S_{pA}$ $(p>1)$ are convex functions of the subsystem size, only linear at $\beta=0$, more convex at larger $|\beta|$ and larger than those of thermal mixed states at the same energy density. This means that $\textrm{Tr}_{A}{\left(\hat{\rho}_{A}\right)^{p}}$ is exponentially small with the subsystem size. The states we are dealing with here are mixed states, and generally $S_{pA}\neq S_{pB}$, unless the partition is in the middle. For a thermal state at temperature $\beta^{-1}$, the Rényi entropy for subsystem $A$ contains contributions from both entanglement and thermal mixture. Refs. Korepin (2004); Calabrese and Cardy (2004, 2009) discuss the crossover of Rényi entropy from the logarithmic scaling ($l_{A}/\beta\ll 1$) to the volume law ($l_{A}/\beta\gg 1$) using finite temperature CFT. This finite temperature behavior can place an additional challenge in the preparation of ground states of cold atom systems in optical lattices Unmuth-Yockey et al. (2017). If $A$ equals the whole system, $S_{pA}$ is the diagonal (thermal) Rényi entropy for the DE (CE). In particular, for the time-evolved density matrix $\hat{\rho}^{t}$, taking the second order Rényi entropy as an example, it can be expressed by

[TABLE]

where $\sum_{\{m,n,i\}}$ is a sum over $m$, $m^{\prime}$, $n$, $n^{\prime}$, $i_{B}$, and $i^{\prime}_{B}$, and where $i_{B}$ and $i_{B}^{\prime}$ run over all the states of a basis for subsystem $B$. If $A$ is the whole system (size $L$), then $B$ has size zero, and $n=m^{\prime}$ and $n^{\prime}=m$, and the phase is therefore always zero. The same argument holds for all orders of the Rényi entropy. The Rényi entropy for the whole system is therefore conserved during the time evolution.

III Weights in Ensembles

In global quench protocols, the energy uncertainty $\Delta E/\bar{E}$, which only depends on the DE, is algebraically small with the system size ($1/\sqrt{L}$) for pure initial states Rigol et al. (2008), which is also numerically verified for our case. To show this explicitly, writing the Hamiltonian after quench as $\hat{H}_{F}=\hat{H}_{0}-\hat{H}_{W}$ and following the derivation in Ref. Rigol et al. (2008), we calculate the energy width of the diagonal ensemble, which is equal to the variance of energy in the initial thermal state,

[TABLE]

If $\hat{H}_{W}$ is a sum of local operators with finite matrix elements, and in the absence of long-range correlations, $\Delta E\sim\sqrt{L}$ still holds using the same argument in Rigol et al. (2008). In particular, in our case, $\hat{H}_{W}$ is a large chemical potential on the right half of the lattice, and the particles are confined in the left half in the initial state, so the last three traces are zero in Eq. (14). So $\Delta E$ equals the variance of the initial energy before quench in a thermal state, which obviously scales as $\sqrt{L}$. The small energy window $[\bar{E}-\Delta E,\bar{E}+\Delta E]$ defines a microcanonical ensemble which is equivalent to the CE in the thermodynamic limit. As our initial states already have the Boltzmann weights, it is interesting to explore and compare the weights in the whole spectrum in the DE, $W^{d}_{n}$, and those in the corresponding CE, $W^{c}_{n}$, and see if these spectra will match by changing the initial temperature. Similar comparisons have been conducted in studies of global quantum quenches from pure states Sorg et al. (2014) and thermal states Zhang et al. (2011); He and Rigol (2012).

In Fig. 1, we show the 2-d histograms for natural logarithms of these weights in the DE [Figs. 1(a), 1(c) and 1(e)] and the CE [Figs. 1(b), 1(d) and 1(f)]. We choose $100\times 100$ bins, and the color scale on each bin represents the number of states inside it. The results are for spinless fermions expanding from thermal states with three different initial temperatures. We see that in the DE, the weights lie in a narrow band for all temperatures. There is a tail bending down for the spectrum at low energies, but with a much smaller density of states. The weights in the CE are all straight lines as should be the case by definition. We have included these as dashed lines in the left column of Fig. 1 to compare the slopes. At low initial temperature, e.g. $\beta_{0}=1$, the slopes of the narrow straight band in the DE and the line for the CE are very different. As we increase the initial temperature, they match better and better. It can be understood in the following way. Note that the summation in Eq. (8), $\sum_{m}E_{m}W^{d}_{m}$, can be divided into two parts: one is the summation in the bending tail of the weights in the DE which is close to the bottom of the energy spectrum, the other one is the summation in the remaining narrow straight band. Although the tail has small density of states, it can still dominate the summation at very low initial temperatures, because the weight in the DE, $W^{d}_{m}=\sum_{i}\rho^{0}_{i}\langle m|i\rangle\langle i|m\rangle$, contains the Boltzmann weight $\rho^{0}_{i}=e^{-\beta_{0}e_{i}}/Z_{0}(\beta_{0})$ given by the initial temperature. As the tail is bending down, the line of the CE must rotate anticlockwise from the straight narrow band of the DE to fit the data in it at low initial temperatures [see Fig. 1(e)]. For higher initial temperatures, the narrow straight band dominates the summation, and Eq. (8) is effectively a linear fit for this band. Note that for the highest initial temperature $\beta_{0}=0.1$, the effective temperature is negative, and the weights in the DE and the CE match perfectly in the densest part of the spectrum [see Fig. 1(a)]. This match between the DE and the CE weights results in agreements for expectations of almost all equithermal operators Garrison and Grover (2018) in the two ensembles and the model has strong thermalization Bañuls et al. (2011) in the negative temperature regime.

The results for bosons are depicted in Fig. 2. The weights in the DE also lie in a straight narrow band. The tail in the low energy spectrum only bends down slightly. So the difference in slopes for the DE and the CE is much smaller than in the case of fermions. We still see that the higher the temperature, the better is the match between the DE and the CE. The slopes also agree accurately at very high initial temperature $\beta_{0}=0.01$. Note that a much higher initial temperature is needed for the effective temperature to be negative. We will see that the negative effective temperature actually does not exist in the thermodynamic limit for bosons, but it does for fermions. This is discussed in Section VII.

IV Entropy For Ensembles

In this section we investigate the entropy profiles for the initial state, the state after long-time expansion, the DE density matrix and the CE density matrix. We calculate the quantity in Eq. (12) for $n=1$ (the von Neumann entropy) and $n=2$ (the second order Rényi entropy), and subsystems containing left $l_{A}=1,2,...,20$ sites. We choose the time evolved density matrix to be at a long time after expansion, $tJ=6000$, where thermalization has occured for a long time. Numeric results show the equilibration of entropy after $tJ=50$ for all cases considered here. We compare the results for different temperatures for both spinless fermions and bosons.

IV.1 Von Neumann Entropy

Fig. 3 shows the von Neumann entropy of the reduced density matrix for the subsystem containing left $l_{A}$ sites. The results are for spinless fermions expanding from thermal states with four different initial temperatures. For initial thermal states, the particles are initially confined in the left half part, so the Rényi entropy is identical for all subsystems bigger than half. We can clearly see from Fig. 3(a) the volume law in $S_{1}$ for up to six subsystems at high temperature $\beta_{0}=0.1$. In Fig. 3(b) a lower temperature $\beta_{0}=0.4$ results in five subsystems satisfying the volume law. As we decrease the initial temperature further, there are only three subsystems satisfying volume law in Fig. 3(c) for $\beta_{0}=1$ and no sign of volume law in Fig. 3(d) for $\beta_{0}=10$. These results show the crossover from the volume law for subsystems at high temperatures to the logarithmic scaling in the ground state. The CE density matrices all show volume law for many subsystems. Because the effective temperature is always high even for very low initial temperature (see Fig. 3(d)). Still, the higher the effective temperature is, the more subsystems satisfying the volume law. The von Neumann entropy increases linearly up to 11 sites in Fig. 3(a) where the negative effective temperature should be considered to be “larger” than the infinite temperature. And it increases linearly up to 6 sites in Fig. 3(d). Intuitively, as the Rényi entropy has contributions from both entanglement and thermal mixture, the former increases at first with increasing subsystem size, but decreases for subsystem sizes bigger than half of the whole system. So it is reasonable that the entropy bends down from the linearly increasing line. For a detailed study of properties of thermal Rényi entropy, see Bonnes et al. (2013).

In Fig. 1 we have shown that the weights in the DE are very similar to the weights in the CE for high initial temperatures. So the diagonal entropy and thermal entropy should be very close. This is confirmed in Fig. 3(a), (b) and Fig. 6(a), (b), where the Rényi entropy of reduced density matrices for $\hat{\rho}_{d}$ and $\hat{\rho}_{c}$ agrees for almost all subsystems. At lower temperatures, on one hand, the deviation for the full systems is consistent with the results in Fig. 1(e)(f) where it shows that the distribution of weights is very different for the DE and the CE. On the other hand, $S_{1A}$ still has no visible difference for subsystems bigger than half even if the distribution of weights is very different for the whole system. This can be seen in the inset of Fig. 5, where the maximal size of the subsystem in which the von Neumann entropy difference between the DE and the CE, $\Delta S_{1cd}=|S_{1}^{c}-S_{1}^{d}|$, is smaller than $0.05$ as a function of $\beta_{0}$ is plotted. For $0.2<\beta_{0}<1.25$, roughly we have $N_{\Delta S_{1}<0.05}\sim\log_{10}{T_{0}}$. In the main plot of Fig. 5, we plot $N_{\Delta S_{1}<0.003}$ as a function of $\beta_{0}$, where $N_{\Delta S_{1}<0.003}\sim\log_{10}{T_{0}}$ for $0.1<\beta_{0}<0.7$. Note that it can be proved that $S_{1}$ of the DE for the whole system only has non-extensive difference from that of the microcanonical ensemble, no matter whether the initial state is pure or mixed D’Alessio et al. (2016). And this non-extensive difference can be decreased by increasing the initial temperature in our case.

Finally, we check the reduced von Neumann entropy for the time-evolved density matrix in Fig. 3. We see that at the high temperature $\beta_{0}=0.1$, the time-evolved density matrix has the same $S_{1}$ as the CE and the DE density matrices for subsystems up to eleven sites. And ten-, eight- and six-site subsystems have no visible difference in von Neumann entropy for temperatures $\beta_{0}=0.4$, $\beta_{0}=1$ and $\beta_{0}=10$, respectively. The maximal size of the subsystem in which the difference between $S_{1}$ for $\hat{\rho}_{t}$ at $tJ=6000$ and $S_{1}$ for the CE, $\Delta S_{1tc}=|S_{1}^{t}-S_{1}^{c}|$, is smaller than $0.003$ ($0.05$) as a function of $\beta_{0}$ is plotted in the main (inset) figure of Fig. 5. The results for $\Delta S^{av}_{1tc}=|\bar{S_{1}^{t}}-S_{1}^{c}|$ with average von Neumann entropy $\bar{S_{1}^{t}}$ over $tJ=81,82,...,200$, are also plotted. The plots for $\Delta S_{1tc}$ and $\Delta S^{av}_{1tc}$ are almost the same in the inset figure, indicating that the time fluctuations are small enough that the instantaneous entropy agrees with the time-averaged entropy at the precision of $0.05$. It can also be seen that roughly $N_{\Delta S_{1}<0.05}\sim\log_{10}{T_{0}}$ for $0.4<\beta_{0}<5$. Nevertheless, the plots for $\Delta S_{1tc}$ and $\Delta S^{av}_{1tc}$ in the main figure are different because $0.003$ is smaller than the fluctuations of $S_{1}^{t}$. In this case $N_{\Delta S^{av}_{1tc}<0.003}\sim\log_{10}{T_{0}}$ only for $0.1<\beta_{0}<0.7$. It is concluded that when the initial temperature is higher, the maximum subsystem size that has the same von Neumann entropy in $\hat{\rho}_{t}$ and $\hat{\rho}_{d}$ as that in $\hat{\rho}_{c}$ is bigger and is expected to scale as $\log_{10}{T_{0}}$. And it is always smaller for $\hat{\rho}_{t}$ than for $\hat{\rho}_{d}$ because $\hat{\rho}_{t}$ preserves the initial entropy. Note that the reduced $S_{1}$ increases first to the highest point and then bends down until the full von Neumann entropy goes back to the initial full von Neumann entropy. This is because the full von Neumann entropy is invariant under the unitary change of basis, which is the unitary time evolution here.

The results for bosons depicted in Fig. 4 are similar to those for spinless fermions. They are summarized as follows. For initial thermal state, there are six, six, five and three sites having volume law von Neumann entropy at initial inverse temperature $\beta_{0}=0.01$ [Fig. 4(a)], $\beta_{0}=0.1$ [Fig. 4(b)], $\beta_{0}=1$ [Fig. 4(c)] and $\beta_{0}=10$ [Fig. 4(d)] respectively. And for the canonical ensemble density matrix, the von Neumann entropy in the left eleven, eleven, ten and seven sites show volume law at these four initial temperatures. The DE von Neumann entropy is identical to the CE von Neumann entropy for all subsystems at initial inverse temperatures $\beta_{0}=0.01,0.1$. They are equal in the left fifteen sites and left eight sites at initial inverse temperatures $\beta_{0}=1$ and $\beta_{0}=10$, respectively. These results are also consistent with those in Fig. 2. The time-evolved density matrix at $tJ=6000$ has the same von Neumann entropy with both $\hat{\rho}^{d}$ and $\hat{\rho}^{c}$ in the left eleven, eleven, nine and six sites in the four graphs respectively.

IV.2 The second order Rényi entropy

Since the results for the two models are very similar, we only consider spinless fermions in this section. As shown in Fig. 6, the second order Rényi entropy also has high agreement for the CE and the DE, that is, nineteen, eighteen, eleven and eight sites have the same value for $\beta_{0}=0.1,0.4,1,10$, respectively. And the time-evolved state has twelve, eleven, ten and eight sites with $S^{t}_{2A}$ equal to them. The main difference from Fig. 3 and Fig. 4 is that, unlike the von Neumann entropy where the higher the initial temperature is, the bigger the subsystems having linearly increasing entropy are, there are five, six, twelve (sixteen for the DE) and six sites having linearly increasing $S^{d(c)}_{2A}$, respectively. We find that the entropy bends up for high initial temperatures [the data of the DE and the CE lie on or above the black dotted linear fit in Fig. 6(a)]. This convexity is weakened by decreasing the initial temperature and finally the entropy becomes concave again [the data of the DE and the CE lie on or below the black dotted linear fit in Fig. 6(d)]. As mentioned in Ref. Lu and Grover (2019), $S_{pA}$ $(p>1)$ are convex functions of the subsystem size, which is verified in our models with twenty sites and quarter filling, but just in a small energy window around the infinite temperature eigenstate. Moreover, more obvious convexity is observed in high energy levels near the top of the spectrum while the levels near the bottom of the spectrum are all concave, which explains why $S^{c(d)}_{2A}$ with negative effective temperature is more convex in Fig. 6 even if $|\beta_{eff}|$ is smaller. Note that at low temperatures, the second order Rényi entropy of the DE can be bigger than that of the CE for big subsystems, while the von Neumann entropy of the DE is always no greater than that of the CE.

As the second order Rényi entropy can be measured experimentally by a beam-splitter operation implemented via a controlled tunneling operation between the two copies of a many-body state Daley et al. (2012); Islam et al. (2015), we investigate the details of time evolution of $S^{t}_{2A}$ in this section. Note that for twenty sites and quarter filling, the numerical results show small time fluctuations in Rényi entropy for all subsystems at long times, even for the lowest temperature considered here $\beta_{0}=10$, which indicates that the time-dependent terms are cancelled out due to dephasing. In Eq. (13), the only time-independent terms are those with $m=n,m^{\prime}=n^{\prime}$ or $m=n^{\prime},m^{\prime}=n$. So for long times,

[TABLE]

where $T_{AB}^{mm^{\prime}}=\textrm{Tr}_{A}{\left[(\textrm{Tr}_{B}{|m\rangle\langle m|})(\textrm{Tr}_{B}{|m^{\prime}\rangle\langle m^{\prime}|})\right]}$ and $T_{BA}^{m^{\prime}m}=\textrm{Tr}_{B}{\left[(\textrm{Tr}_{A}{|m^{\prime}\rangle\langle m^{\prime}|})(\textrm{Tr}_{A}{|m\rangle\langle m|})\right]}$. The second summation contributes to the difference between the equilibrated $S_{2A}$ and the $S_{2A}$ for the DE. If it is negligible compared to the first summation, the time-evolved Rényi entropy equilibrates to the value for the DE which agrees with the value for the CE. For small subsystem sizes $l_{A}\ll l_{B}$, consider any $m\neq m^{\prime}$, we expect $T_{AB}^{mm^{\prime}}\gg T_{BA}^{m^{\prime}m}$ and $S^{t}_{2A}\approx S^{d}_{2A}$ because the overlap of $\textrm{Tr}_{B}(|m\rangle\langle m|)$ for two $m$ values is much larger than the corresponding overlap for $\textrm{Tr}_{A}(|m\rangle\langle m|)$. Furthermore, as $\sum_{m,m^{\prime}}\rho^{0}_{mm}\rho^{0}_{m^{\prime}m^{\prime}}=1$, the first summation is equal to some typical value of $T_{AB}^{mm^{\prime}}$, denoted as $T_{AB}^{o}$, while the second summation can be written as

[TABLE]

where $T_{BA}^{o}$ is some typical value of $T_{BA}^{m^{\prime}m}$, $S_{2}^{\beta_{0}}$ is the second order Rényi entropy for the initial state, $S_{2}^{d}$ is the second order Rényi entropy for the DE. $S_{2}^{d}$ is extensive with $L$ because it is close to $S_{2}^{c}$ (see Fig. 6), while $S_{2}^{\beta_{0}}$ is extensive with $L/2$. Because the average entropy per site for quarter filling is bigger than that for half filling and the effective temperature is higher than the initial temperature after the expansion, the second term in Eq. (16) is negligible for large $L$. So as long as the initial state is thermal, Eq. (16) is at least exponentially small with $L/2$.

Note that the dependence of $T_{AB}^{mm^{\prime}}$ and $T_{BA}^{m^{\prime}m}$ on the initial temperature and the system size also contributes to the difference between the time-evolved Rényi entropy and the DE Rényi entropy. If Eq. (16) is small,

[TABLE]

where $T^{o}_{BA}<\textrm{Tr}_{B}\left[\left(\textrm{Tr}_{A}|m_{o}\rangle\langle m_{o}|\right)^{2}\right]$ is also exponentially small with the subsystem size for $|m_{o}\rangle$ being a finite-energy density eigenstate of a chaotic system Lu and Grover (2019). We plot Figs. 7 and 8 to quantitatively investigate the second term in Eq. (17). To plot Fig. 7, we calculate the time evolution of $S_{2A}^{t}$ for subsystems containing the left $9,10,11$ sites respectively to $tJ=200$, with time step size $\tau J=1$. To avoid large deviations due to the fluctuations in time and the short-time far-from-equilibrium states, we average $S_{2A}^{t}$ using only the last $120$ steps of the time evolution data. In Fig. 7, we plot the difference between $S_{2A}^{d}$ for the DE and the time-averaged value $\Delta S_{2A}=S_{2A}^{d}-\bar{S}_{2A}$ as a function of the second order Rényi entropy for the initial thermal state $S^{\beta_{0}}_{2}$. Varying the number of long-time data points in the time-average calculation does not change the plot. The exponential decay of $\Delta S_{2A}$ with the initial thermal entropy can be clearly seen from these numerical results. The decay rates are almost the same for three subsystems at least for $\beta_{0}>0.2$. The results at $tJ=6000$ are also plotted, which are very close to the results for time-averaged entropy. The inset of Fig. 7 shows $S^{\beta_{0}}_{2}$ as a function of initial temperature $T_{0}$, where we see $S^{\beta_{0}}_{2}$ is firstly linear with $T_{0}$ at low initial temperature and saturates at high initial temperature. So $\Delta S_{2A}$ also exponentially decays with the initial temperature when $T_{0}<1$. In Fig. 8, the same quantity for $l_{A}/L=1/2$ is plotted as a function of the system size, where we see that $\Delta S_{2A}$ is exponentially small with the system size. The decay rate is higher for higher initial temperatures. The inset (a) confirms that $S_{2A}^{\beta_{0}}$ is proportional to the size of the system. So $\Delta S_{2A}$ also exponentially decays with $S_{2A}^{\beta_{0}}$ across different system sizes. The inset (b) shows the linearity of $S^{d}_{2}$ with the system size, where the weak convexity is due to the increase of effective temperature for larger system sizes. The two insets also confirm that $e^{-S_{2}^{d}}$ is negligible in Eq. 16, at least for $L\geq 12$ where the data of $S^{d}_{2}$ is more linear with $L$. The results at $tJ=6000$ are not close to the time-averaged ones due to the large time fluctuations for small systems.

Finally, we investigate the time fluctuations of the Rényi entropy. We calculate the variance of $S^{t}_{2A}$ in time,

[TABLE]

where we take $N_{\tau}=120$ and $\{\tau_{i}\}=\{81,82,...,200\}$. The variance as a function of $S^{\beta_{0}}_{2}$ is depicted in Fig. 9, where we see that the fluctuations also exponentially decay with $S^{\beta_{0}}_{2}$, at least for not very high temperatures ($\beta_{0}=0.1$ for the largest $S^{\beta_{0}}_{2}$ in the figure). The inset also shows that the variance is exponentially small with the system size, with a bigger decay rate for higher temperatures.

V Reduced Density Matrices

The above observations indicate that the reduced density matrices of $\hat{\rho}^{t}$ at long time, $\hat{\rho}^{d}$ and $\hat{\rho}^{c}$ may be very similar, for not too large subsystems. In order to check this, we calculate the eigenvalues of the reduced density matrices $\lambda_{A}$ and plot $\ln{\lambda_{A}}$ as a function of the eigenenergy $E^{m}_{A}$ of the Hamiltonian associated with the subsystem $\hat{H}_{A}$. Another important question is to decide if these eigenvalues $\lambda_{A}$ are equal to the weights in thermal density matrix of the subsystem $\hat{\rho}^{c}(A)=e^{-\beta_{eff}\hat{H}_{A}}/Z_{A}$. In order to answer this question, we consider subsystems containing more than five sites, for which the reduced density matrix can be divided into six sub-sectors which have zero, one, two, three, four and five particles respectively. If the sub-sector contains $p$ particles, the bath of the subsystem should have $5-p$ particles. The reduced density matrices of density matrices considered above are all mixtures of states with fixed number of particles in subsystem $A$. We need to construct the thermal state of subsystem $A$ with a correct mixture of thermal states in these six sub-sectors. Notice that the Hamiltonian can be written as $\hat{H}=\hat{H}_{A}+\hat{H}_{B}+\hat{H}_{AB}$. $\hat{H}_{A}$ and $\hat{H}_{B}$ are extensive quantities, while $\hat{H}_{AB}$ is a non-extensive local interaction with negligible contribution to the total energy. Then the product of each pair of eigenstates $|E_{A}\rangle,|E_{B}\rangle$ can approximate an eigenstate of $H$ with energy $E\simeq E_{A}+E_{B}$. A thermal state can then be approximated by a product of two thermal states associated with the subsystem and its complement. To be specific, let $\{|E_{A}^{m_{p}}\rangle\}$ be the eigenstates of $H_{A}^{p}$ with $p$ particles, $\{|E_{B}^{m_{5-p}}\rangle\}$ be the eigenstates of $H_{B}^{5-p}$ with $5-p$ particles, the whole thermal density matrix of two uncorrelated subsystems $A$ and $B$ at temperature $\beta_{eff}$ is

[TABLE]

where $Z_{AB}=\sum_{p=0}^{5}Z_{A}^{p}Z_{B}^{5-p}$ is the partition function of the whole system. So the thermal density matrix of $A$ is

[TABLE]

Thus we have obtained the thermal state of subsystem $A$ as a weighted mixture of states with fixed number of particles. The weights can be obtained by computing the spectrum in subsystems $A$ and $B$ separately.

Fig. 10 shows the results for the subsystem containing left six sites. Fig. 10(a) is for spinless fermions with $\beta_{0}=0.1$ and Fig. 10(b) is for bosons with $\beta_{0}=0.01$. Fig. 11 depicts the results for spinless fermions with $\beta_{0}=0.1$ [Fig. 11(a)] and $\beta_{0}=0.4$ [Fig. 11(b)], where the subsystem contains left ten sites. And Fig. 12 shows the results for bosons with $\beta_{0}=0.01$ [Fig. 12(a)] and $\beta_{0}=1$ [Fig. 12(b)], where the subsystem also contains left ten sites. Fig. 10, Fig. 11(a) and Fig. 12(a) have negative effective temperatures, so the slopes of the lines are positive, while the effective temperatures are positive for Fig. 11(b) and Fig. 12(b), the slopes are negative. The points linearly arranged belong to the same sub-sector with a certain particle number. From top to bottom, they are zero, one, two, three, four, five particle states, respectively. First note that for spinless fermions the spectrum of reduced CE, $\textrm{Tr}_{B}{(\hat{\rho}^{c})}$, match perfectly with that of $\hat{\rho}^{c}(A)$ for all cases shown here. However, there are two main deviations for bosons. First, there are some ripples in the spectra of reduced CE around that of $\hat{\rho}^{c}(A)$ for sub-sectors containing more than two particles. This effect is more obvious for the lower temperature $\beta_{0}=1$ than for the higher temperature $\beta_{0}=0.01$ in the same subsystem, see Fig. 12. Second, there are discontinuities near the top of the spectrum for sub-sectors containing more than two particles. The states near the top of the spectrum in those sub-sectors jump higher from the lines for the negative effective temperature but jump lower from the lines for the case of positive effective temperature. But overall they still match well.

For spinless fermions, it can be seen that in Fig. 10(a), the four density matrices have almost the same spectra, except that the line of the five-particle sector in the time-evolved state rotates anti-clockwise a little. In Fig. 11(a), the initial temperature is the same, but the subsystem is bigger, then only the first two sub-sectors containing zero and one particle have perfect match for all four density matrices. While for the sub-sector containing two particles, the eigenvalues of $\textrm{Tr}_{B}{\hat{\rho}^{d}}$ and $\textrm{Tr}_{B}{\hat{\rho}^{t}}$ become smaller near the bottom of the spectra, similar to Fig. 1 where the DE also has lower weights near the bottom of the spectra. For sub-sectors with more particles, the reduced DE can always stay around the reduced CE because the similarity of weights between DE and CE as shown in Fig. 1, while $\textrm{Tr}_{B}{\hat{\rho}^{t}}$ has bigger and bigger deviations. In Fig. 11(b) where the initial temperature is lower, the spectra in the first two sectors still match well. For the sub-sector with two particles, the eigenvalues for $\textrm{Tr}_{B}{\hat{\rho}^{d}}$ and $\textrm{Tr}_{B}{\hat{\rho}^{t}}$ are smaller near both the top and the bottom of the spectrum, also similar to the behavior of weights in Fig. 1. The fact that the reduced DE have bigger deviations at lower initial temperature is also consistent with the result in Fig. 1 that the DE has bigger difference from the CE at lower initial temperature. The spectrum of $\textrm{Tr}_{B}{\hat{\rho}^{t}}$ is also very different from the other three for the five-particle states. Note that the big deviations in the time evolved state is not due to time fluctuations. We have checked the results for other ten different times and they behave almost the same. So the eigenvalues of $\textrm{Tr}_{B}{\hat{\rho}^{t}}$ still equilibrate but not to those of reduced DE.

For bosons, as shown in Fig. 12, we can still see the ripples and discontinuities in $\textrm{Tr}_{B}{(\hat{\rho}^{d})}$ and $\textrm{Tr}_{B}{(\hat{\rho}^{t})}$ for the sub-sectors containing more than two particles. The spectrum for reduced DE match well to the reduced CE for all three cases considered here, because the DE and the CE for bosons are very similar for initial inverse temperature smaller than $1$, as shown in Fig. 2. For $\beta_{0}=1$, Fig. 2(c) shows that most of the weights in DE are a little smaller than those in CE, which also happens here for the reduced DE and the reduced CE. $\textrm{Tr}_{B}{(\hat{\rho}^{t})}$ have eigenvalues matching well to the other three for sub-sectors containing no more than four particles, but big differences appear for five-particle states in the subsystem with ten sites.

For all the cases considered here, most of the eigenvalues of all four density matrices in at least the first four sub-sectors with no more than three particles agree with each other. The contribution from these eigenvalues dominate the value of Rényi entropy, so the entropy of the subsystem containing the left ten sites or fewer all have indiscernible difference in Sec. IV. Another question is if the corresponding eigenstates, $|\lambda_{A}\rangle$, are the same with the eigenstates of the Hamiltonian associated with the subsystem, $|E_{A}^{m}\rangle$. If it is, then any observables residing in this subsystem are expected to thermalize as long as they depend little on the high energy eigenstates belonging to the sub-sectors with large number of particles. To check this, we plot the overlaps between eigenstates of the three density matrices and those of $\hat{H}_{A}$ in Fig. 13. The results are for spinless fermions with $\beta_{0}=0.1$ [Fig. 13(a-f)] and $\beta_{0}=0.4$ [Fig. 13(g-l)]. Fig. 13(a-c)(g-i) are for the subsystem containing the left six sites and Fig. 13(d-f)(j-l) are for the subsystem containing the left ten sites. Fig. 13(a,d,g,j) are for the reduced time-evolved states, Fig. 13(b,e,h,k) for the reduced DE and Fig. 13(c,f,i,l) for the reduced CE. For $\beta_{0}=0.1$, we see from Fig. 13(c)(f) that most of the eigenstates of the reduced CE are the same with those of $\hat{H}_{A}$ (referred as “good” eigenstates here). And we have shown that the eigenvalues are very close. So the reduced density matrix of a thermal state is still thermal with the same temperature. This is the intensive property of temperature in quantum mechanics. The reduced DE for the subsystem containing six sites still have many good eigenstates, most of which are in the sub-sectors with no more than two particles, while for the subsystem containing ten sites it has only a small portion of eigenstates that are good. As the expectation value of an observable depend on the diagonal ensemble, we conclude that nearly all obervables thermalize as long as they reside in small subsystems and few-particle sub-sectors. The reduced time-evolved state has much fewer good eigenstates. But it can be seen that for the subsystem containing six sites, the first two sub-sectors containing no more than one particle (first 7 states) still have good resutls. So for small subsystems and few-particle sub-sectors, the time-evolved states becomes thermal at long enough time. The results for the lower initial temperature $\beta_{0}=0.4$ are better than those for $\beta_{0}=0.1$, where it is clear that there are more good eigenstates in both the reduced time-evolved state and the reduced DE. Note that for the six-site subsystem, unlike the case with $\beta_{0}=0.1$ where most of the good eigenstates are in low particle number sub-sectors, $\beta_{0}=0.4$ has most of the good eigenstates in high particle number sub-sectors, which may due to the different signs of effective temeprature. We also investigate the results for $\beta_{0}=1$, which have similar behaviors to $\beta_{0}=0.4$ but are a little worse than $\beta_{0}=0.4$. Note that more good eigenstates does not mean that the density matrix is closer to thermal. As high particle number eigenstates have much lower weights, the reduced time-evolved states for higher initial temperatures may be still closer to the thermal states.

We conclude this section with a comment on the theory of equilibration and thermalization. Intuitively, few-particle sub-sectors have more degrees of freedom serving as their bath, so they should behave more like the corresponding thermal states. For the sub-sector containing five particles, the bath only has one degree of freedom. But there still are a large number of dephasing terms, which could result in the equilibration of eigenvalues of reduced time-evolved states but not to those of reduced DE, similar to the behaviors of Rényi entropy. In Ref. Linden et al. (2009), a theorem states that the time averaged distance between the reduced time-evolved density matrix and the reduced DE density matrix is bounded by $1/2\sqrt{d_{A}/d^{eff}(\textrm{Tr}_{A}{(\hat{\rho}_{d})})}\leq 1/2\sqrt{d_{A}^{2}/d^{eff}(\hat{\rho}_{d})}$, where $d_{A}$ is the dimension of the subsystem, $d^{eff}(\textrm{Tr}_{A}{(\hat{\rho}_{d})})$ is the effective dimension of the reduced DE and $d^{eff}(\hat{\rho}_{d})$ is the effective dimension of the DE. From the stronger bound we see that the subsystem must be smaller than half in the sense $l_{A}/L<1/2$. For the weaker bound, the maximal $d^{eff}(\hat{\rho}_{d})$ is the dimension of the Hilbert space for the whole system when the weights in the DE are constant. In our case, the weights in the CE and the DE are very similar at high effective temperatures, so for $l_{A}/L<1/2$, the subsystem equilibrates to the reduced DE density matrix in the thermodynamic limit. Ref. Linden et al. (2009) also proves that the effective dimension of the diagonal ensemble is the order of the dimension of the Hilbert space (with exponentially small probability to be smaller). We can directly count the dimension of the whole system and the subsystem, and use Stirling formula to estimate the bound. When $l_{A}/L=1/2$ and large $L$, the dimension of the whole system for spinless fermions $\ln{d^{f}_{w}}\approx(2\ln{2}-3/4\ln{3})L$, and $\ln{d^{f}_{A}}\approx L/2\ln{2}$. So $(d^{f}_{A})^{2}/d^{f}_{w}\approx e^{0.131L}$ is exponentially large with the size of the system. A careful calculation shows that we need $l_{A}/L<0.4175$ for the weaker bound to be exponentially small at high temperature. Similarly for bosons, $\ln{d^{b}_{w}}\approx(5/4\ln{5}-2\ln{2})L$, $\ln{d^{b}_{A}}\approx(3/4\ln{3}-1/2\ln{2})L$, $(d^{b}_{A})^{2}/d^{b}_{w}\approx e^{0.33L}$. And we need $l_{A}/L<0.2045$ for the weaker bound to be exponentially small. The numeric results for $l_{A}/L=1/2$ show good agreement between eigenvalues but not good agreement between eigenstates. Following the derivation in the Appendix A of Ref. Linden et al. (2009), replacing the coefficients by the matrix elements considered here, the average distance is bounded by $1/2\sqrt{d_{A}e^{-S_{2}^{\beta_{0}}}T_{BA}^{o}}$, where again $T_{BA}^{o}$ is a typical value of $T_{BA}^{mn}$. For $\beta_{0}=0$, $d_{A}e^{-S_{2}^{\beta_{0}}}$ is at most an algebraically increasing function of system size, while $T_{BA}^{o}$ is typically exponentially small. So in our case, it is possible for subsystems bigger than half of the system ($l_{A}/L>1/2$) to equilibrate to the corresponding reduced DE density matrix in the thermodynamic limit.

Next, the DE is a weighted sum of eigenstate projectors. ETH in the strong sense states that every eigenstate is thermal for few-body observables Kim et al. (2014). And the similarity of reduced density matrices between the DE and the CE can be explained by the strong form of ETH which states that the reduced density matrix of a single finite energy density eigenstate of chaotic many-body quantum systems is equivalent to the thermal density matrix of the subsystem as long as the subsystem is much smaller than its complement and as long as $l_{A}/L<1/2$ for many observables Garrison and Grover (2018). The subsystem ETH Dymarsky et al. (2018) states that the norm of the “off-diagonal” matrix $\hat{\rho}_{A}^{mn}=\textrm{Tr}_{B}{(|m\rangle\langle n|)}$ ($m\neq n$) is exponentially small for subsystems with $l_{A}/L<1/2$, which explains the similarity between reduced $\hat{\rho}^{t}$ and reduced DE. Inside the small window where the energy density is peaked, all the reduced density matrix of eigenstates give the same thermal states for the subsystem. We show that even the full density matrix of the DE is close to the CE at high temperatures for the models we study here. So at high temperatures, the size of the subsystem that thermalizes just depends on the size of the subsystem that equilibrates to the reduced DE.

VI Momentum Distribution Function

In this section, we use the momentum distribution function (MDF) to test the above conclusions. The MDF for fermions is

[TABLE]

where $p,q=1,2,...,L$ are the real-space locations. And the MDF for bosons is

[TABLE]

The MDF is a few-body observable because it contains a summation of two-body observables, so it is expected to thermalize according to the above results. Regardless, we calculate the MDF for the initial state $\hat{\rho}^{0}$, the time-evolved state $\hat{\rho}^{t}$ at long time $tJ=6000$, the DE density matrix $\hat{\rho}^{d}$ and the CE density matrix $\hat{\rho}^{c}$. If the MDF equilibrates, its expectation value using $\hat{\rho}^{t}$ should be the same as that using $\hat{\rho}^{d}$. And if the MDF thermalizes, all the expectation values, using $\hat{\rho}^{t}$, $\hat{\rho}^{d}$ and $\hat{\rho}^{c}$ should be the same.

Fig. 14 shows the results for spinless fermions with the same four initial temperatures as those in Fig. 3. Let’s first briefly describe the MDF in the initial state. At low initial temperature $\beta_{0}=10$ [Fig. 14(d)], there are plateaus around $k=0,\pm 1.8$. The momenta around $k=0,\pm\pi$ have the highest population, while the momenta around $k=\pm 1.8$ have the lowest population. These behaviors can be understood by Fourier transforming the kinetic terms in Hamiltonian Eq. (2), $-2(\cos{k}+\cos{2k})$, which has one global minima $k=0$, two local minima $k=\pm\pi$ and two global maxima $k\approx\pm 1.8$. At low temperature, the global minimum $k=0$ is occupied first, then the nearby levels, one particle per level, due to the Pauli exclusion principle. The height of the plateau is $0.5$ because of normalization by the full volume. If the surrounding levels have higher energy than the other two local minima $k=\pm\pi$ do, particles start to occupy $k=\pm\pi$ and their surrounding levels. Very few particles will occupy the two maxima $k=\pm 1.8$ and their surroundings. At higher temperatures (Fig. 14(a)(b)(c)), the plateaus are smoothed out by thermal fluctuations. Close to the infinite temperature limit, more and more particles occupy the high energy levels, and eventually the MDF becomes flat.

The equilibration and thermalization of the MDF in Fig. 14 are very clear. In Fig. 14(a), $\beta_{0}=0.1$, the MDFs for $\hat{\rho}_{t}$, $\hat{\rho}_{d}$ and $\hat{\rho}_{c}$ do not have a visible difference, indicating strong thermalization of the MDF. The inversion of occupation is very interesting, and is a result of the changing of sign of the temperature. The effective final temperature we calculated for this case is indeed negative. The difference in three MDFs is still extremely small in Fig. 14(b) where the initial temperature is lower $\beta_{0}=0.4$. Further decreasing the initial temperature results in bigger difference between the MDF for $\hat{\rho}_{c}$ and the other two MDFs. This difference is clear in Fig. 14(d), where the MDF for $\hat{\rho}_{c}$ has lower occupation around $k=\pm\pi$, higher occupation around $k=\pm 1.8$, lower occupation for $k\in(-1.1,-0.7)$ and higher occupation for $k\in(-0.7,0)$. These differences are probably finite size effects because they are much smaller than those in smaller systems with the same particle filling. Since these differences are overall relatively small, the MDF in the CE is still a good approximation for the real MDF at long times. It should be emphasized that although there are three slightly asymmetric points for $k\in(0,0.5)$ in Fig. 14(d), which is due to time fluctuations, the MDFs for DE and time-evolved state are the same for all four initial temperatures.

Fig. 15 shows the results for bosons. As the model we consider for bosons only has nearest-neighbor hopping, the hopping energy only have a single minimum at $k=0$ and maxima at $k=\pm\pi$. So for positive effective temperatures [Fig. 15(b)(c)(d)], the occupation is peaked at $k=0$, while for negative effective temperatures [Fig. 15(a)], the occupation is the lowest at $k=0$. Note that the MDF for $\hat{\rho}^{t}$ and $\hat{\rho}^{d}$ has little difference in Fig. 15(a)(b)(c), and the small difference around $k=0$ in Fig. 15(d) comes from time fluctuations. While they have obvious deviations from that of the CE for all cases, we can see that the MDF of the CE is more peaked at $k=0$ for positive effective temperatures in Fig. 15(c)(d), while its value around $k=0$ is smaller than that of the DE for negative effective temperature in Fig. 15(a). In Fig. 15(b), the effective temperature is close to infinity, MDF of the CE is flatter than that of the DE. Again these deviations are much smaller than those in smaller systems, so we believe they are finite-size effects.

From the above results we conclude that equilibration is much easier than thermalization. The intuitive reason is that the time-dependent terms are expected to cancel out due to dephasing. And it can also be partially understood from the similarity of $\textrm{Tr}_{B}{(\hat{\rho}^{d})}$ and $\textrm{Tr}_{B}{(\hat{\rho}^{t})}$ discussed in the last section. To see it more rigorously, we follow the derivation in Ref. D’Alessio et al. (2016); Zhang et al. (2011) and calculate the variance of an observable in time

[TABLE]

where $M$ is the maximal off-diagonal matrix element of $\hat{O}$, then similar to Eq. 16,

[TABLE]

From the same arguments for Eq. (16) we conclude that as long as the initial state is at a finite temperature, the variance exponentially decays with the volume of the system, as long as the off diagonal elements are not exponentially large as $e^{S_{2}^{\beta_{0}}}$. Moreover, even at $T=0$, ETH asserts that the off diagonal elements of few-body observables are exponentially small with the size of the system. So initial thermal states substantially decrease the fluctuations in time.

VII The Effective Temperature

We have seen that the effective temperature obtained by Eq. 8 can be negative. It is interesting that the negative temperature state can be created simply by sudden expansion to a larger volume. The numerical results presented in the previous sections are for finite system sizes. An important question is whether this negative temperature survives in the thermodynamic limit. To answer this question, we do a finite size analysis of the effective inverse temperature $\beta_{eff}$. In order to go to larger system sizes, we consider a configuration where $L/4$ particles expand inside a ring of $L$ sites and carry out exact diagonalization of the final ($L$-site) Hamiltonian. Making use of the translation invariance Sandvik (2010), we can then calculate the spectrum of the final Hamiltonian up to 28 sites with 7 particles for spinless fermions (the dimension of Hilbert space is $1,184,040$) and 24 sites with 6 particles for bosons (the dimension of Hilbert space is $475,020$). We also use finite temperature density matrix renormalization group (DMRG) Verstraete et al. (2004); Feiguin and White (2005); Barthel (2016) with approximated matrix product operators (MPO) applied to matrix product states (MPS) Östlund and Rommer (1995) for larger systems. We perform the calculations using the ITensor C++ Library 111Version 3.0.0, http://itensor.org/, approximating the time evolution operator with MPO to the second order Zaletel et al. (2015) and setting the imaginary time evolution step to be $\tau=0.001$. We calculated systems with $L=60,100,400$ for fermions and $L=28,36,48,64,100$ for bosons. The data for bosons with $L=28,36,48$ are calculated using MPO with periodic boundary conditions, while results of $L=64,100$ are calculated with open boundary conditions. Convergence of finite temperature DMRG results is challenging for large bosonic systems at negative temperatures, and only data that has reliably converged is shown in Fig. 17. For fermionic systems, we use results from exact diagonalization ($L=12,16,20,24,28$) in the curve fitting procedure and the DMRG data for larger systems for checking extrapolation ($L=60,100,400$). These results are shown in Fig. 16. We use all converged data to do curve fitting for bosonic systems. The boundary effects for $L=64,100$ (open boundary conditions) are negligible.

The results for spinless fermions are depicted in Fig. 16. We see, for all initial temperatures shown in the figure, that $\beta_{eff}<\beta_{0}$ and $\beta_{eff}$ decreases as $L$ increases. The effective inverse temperature $\beta_{eff}$ is an approximated linear function of $1/L$. The lines in Fig. 16 are curve fits for data from exact diagonalization ($L=12,16,20,24,28$) with third-degree polynomials. The curve fitting is stable to the degree of the polynomials and has negligible difference when higher degree polynomials are used. The curve fits can be extrapolated and the values of $\beta_{eff}$ in the thermodynamic limit can be extracted. An initial inverse temperature $\beta_{0}=0.2$ already results in a negative temperature state for $L\rightarrow\infty$, and higher initial temperatures result in more negative inverse temperature states after expansion. The validity of the extrapolation with data from exact diagonalization is confirmed by the data from finite temperature DMRG, which are very close and slightly below the extrapolation for $L=60,100$ due to boundary effects (we use open boundary conditions in the DMRG calculations) and coincide with the extrapolation for $L=400$. By repeating the procedure of extrapolation from $\beta_{0}=-0.2$ to $\beta_{0}=0.4$ with stepsize $0.01$, we show the effective inverse temperature as a function of initial inverse temperature in the thermodynamic limit in Fig. 18. It is a smooth function with slight concavity, from where we see that the initial inverse temperature that has zero effective inverse temperature is $0.22$.

Negative temperature occur in systems with bounded spectrum and therefore have a peak in the thermal entropy as function of the energy. For energies above the value where the peak occurs, the thermal entropy decreases with energy and these high energy states have negative temperature. High energy states for the fermionic model correspond to particles occupying neighboring sites (so that the energy is in the form of high interaction energy), or occupying high-energy $k$-states (energy in the form of high kinetic energy), depending on which scale dominates. In either case, and even when the two scales are comparable, having high energy reduces the number of allowed microscopic configurations and the thermal entropy decreases with increasing energy. The thermal entropy is maximum when all the eigenvalues of the density operators are equal, which is the case of infinite temperature in CE. This is the standard picture for thermal entropy and negative temperature. Now, in Sec. IV, we have shown that the three density matrices (exact time-evolved at long times, diagonal, and canonical) all have very similar spectra for small subsystems, and the DE and the CE continue to agree for larger subsystems, even up to the entire system in some cases. So the entanglement entropy of small subsystems behaves the same way as the thermal entropy at negative temperature. Here, we extract the effective temperature by finding a CE with the same total energy as the initial system (Eq. 8), and thus find negative effective temperatures for high energies (high enough initial temperature). In order to draw contrast later to the case of bosons that we discuss below, we note also that even the highest eigenenergy of the fermionic model is extensive and the canonical ensemble is therefore well defined for negative temperature. So as long as the initial energy is higher than the infinite temperature energy of the final Hamiltonian, the effective temperature becomes negative after expansion.

The same plot as Fig. 16 for fermions is shown in Fig. 17 for bosons. As in the case of fermions we see that $\beta_{eff}<\beta_{0}$ but for bosons, $\beta_{eff}$ increases with the system size $L$ for large enough systems. Due to strong finite size effects, we use both exact diagonalization and DMRG data to do curve fitting in order to get an estimate for the error of the extrapolated value of $\beta_{eff}$ in the thermodynamic limit. The solid lines in Fig. 17 are fifth-degree polynomial fits and the estimate for the errors are obtained by changing the degree of the polynomials. For $\beta_{0}=0.1,0.2$, the curve fits are stable to both the degree of polynomials and number of data points included, and the error bars for the extrapolated values at $L\rightarrow\infty$ are too small to be seen in the plot ($<10^{-3}$). This stability lasts down to the extrapolated initial inverse temperature $\beta_{0c}=0.0891(6)$ that has zero effective inverse temperature. The extrapolated $\beta_{eff}$ as a function of $\beta_{0}$ ($\beta_{0c}<\beta_{0}<0.4$) in the thermodynamic limit is depicted in Fig. 18, which exhibits much less concavity compared to the case of fermions. For smaller initial inverse temperatures and negative $\beta_{eff}$, $\beta_{eff}$ decreases at first and then increases for larger systems. The curve fits are very sensitive to the degree of the polynomials used in the curve fitting, so there is large uncertainty in the extrapolated ($L\rightarrow\infty$) value of $\beta_{eff}$ for these cases. Finite temperature DMRG is difficult for large systems of bosons at negative temperature due to high onsite occupation and huge quadratic onsite energy, so we do not have reliable values of $\beta_{eff}$ for $L>48$ with $\beta_{0}=-0.1,-0.2$. We now discuss the upturn in $\beta_{eff}$ with increasing system size $L$ that is exhibited by the bosons but not the fermions. For the bosonic model, since $U/J=3$, the highest energy state corresponds the bosons all occupying a single site, which has energy $\sim N_{p}(N_{p}-1)$, and the next highest energy state has energy $\sim(N_{p}-1)(N_{p}-2)$, where $N_{p}=L/4$ is the total number of bosons. The energy difference between these two states is linear in $N_{p}$. When the temperature is finite and negative and $L$ is large, the Boltzmann weights in CE for the highest energy state dominate Braun et al. (2013), all bosons stay at the same site and the energy is quadratic in $N_{p}$, and therefore quadratic in $L$. Initial states with non-negative or $-1/L$ like initial inverse temperature have extensive total energy, and we expect that the negative temperature should not survive in the thermodynamic limit for bosons with repulsive interactions. However, this needs to be reconciled with the previous argument, based on decreasing entropy with increasing energy, that the final temperature should be negative. The infinite temperature energy of $m$ bosons in $N$ sites is $Um(m-1)/(N+1)$. The energy of the initial Hamiltonian at infinite temperature is $UL(L-4)/8(L+2)$ ($m=L/4,N=L/2$), twice that of the final Hamiltonian at infinite temperature $UL(L-4)/16(L+1)$ ($m=L/4,N=L$). So when the initial state has extensive energy and the initial energy is higher than the energy of the final Hamiltonian at infinite temperature, $\beta_{eff}$ should be negative. These two expectations can be conciliated by having $\beta_{eff}\sim-1/L$ when $\beta_{0}=0$, that is, while $\beta_{eff}$ does indeed become negative for finite systems, it tends back to zero (infinite temperature) in the thermodynamic limit. The form $\beta_{eff}\sim-1/L$ returns the Boltzman weights to their conventional form and ensures that the extensive property of the total energy is maintained. In order to check this last remark, we calculated the total energy for quarter-filled bosonic systems with $L=28,32,36,40,48,56,64,76,88,100$, and inverse temperature given by $\beta=-1/L$ using finite temperature DMRG. The results are plotted in Fig. 19 as a function of system size $L$. The linear fit demonstrates the extensive property of the total energy. For $\beta_{0}=0$, $\beta_{eff}$ is then expected to go to [math] as $\beta_{eff}=-1/L$ for very large systems. Larger scale computations is needed to confirm this numerically, but the statement is consistent with our numerical results within the error bars (Fig. 17). The statement holds for any $0\leq\beta_{0}\leq\beta_{0c}$ in the thermodynamic limit where the initial energy is extensive and exceeds the infinite temperature energy of the final Hamiltonian. Finally, if the initial temperature is already negative, for large $L$ all bosons already form a single particle that moves together, and the total energy is already non-extensive. The effective temperature should then be the same as the initial temperature in the thermodynamic limit. This is consistent with the extrapolation in Fig. 17, within error bars, for $\beta_{0}=-0.2$. In summary, $\beta_{eff}=\beta_{0}$ for $\beta_{0}<0$ and $\beta_{eff}=0$ for $0\leq\beta_{0}\leq\beta_{0c}$, as shown in Fig. 18.

VIII Conclusion

We studied the Joule expansion for one-dimensional nonintegrable quantum lattice systems containing bosons or spinless fermions. We calculated and compared the weights of the eigenstates in the DE and the CE, for different initial temperatures, for the two models. The effective final temperature for the CE is obtained by matching the total energy of the system before and after the expansion. The effective temperature can be negative when the initial temperature is high enough, where the weights of the DE and the CE can match very well. As the DE determines the equilibrated expectation values of observables, it is expected that extensive observables can thermalize at high initial temperatures. We analyzed the von Neumann and second order Rényi entropies in the DE, the CE and for the exact numerical time-evolved state at long times. The agreement between the DE and the CE can extend to the whole system for high initial temperatures because of the similarity in weights. And the agreement between the time-evolved states and the other two can extend beyond $l_{A}=L/2$ at high temperatures. We provided analytical arguments and numerical confirmation that $\Delta S_{2A}$, the difference between the second-order Renyi entropy calculated in the DE and the fully time-evolved time-averaged one at long times, around $l_{A}=L/2$ is exponentially small in system size $L$ (Fig. 8). We examined the eigenvalues and the eigenstates of the reduced density matrices directly and found that the reduced time-evolved density matrix at long times can be thermal for small subsystems. The momentum distribution functions have strong thermalizations at high initial temperature, and show population inversion after the expansion if the effective temperature is negative.

Note again that for initial temperatures larger than a certain value, the system thermalizes in a negative-temperature state. Negative temperature states have been created in cold atoms in spin systemsMedley et al. (2011) and for motional degrees of freedom of a fermionic systemBraun et al. (2013), and are of interest for a number of exotic phenomena that can be realized with themMandt et al. (2011). We propose that negative temperature states can be created by sudden expansion to a larger system size. In this paper we only considered repulsive interactions. For fermions, the negative temperature is a feature that remains in the thermodynamic limit. For bosons with repulsive interactions, there can be no negative temperature state in the thermodynamic limit since the energy spectrum is unbounded. For finite, experimentally realistic system sizes, there can be negative-temperature states and these can also be obtained through Joule expansion.

There are still a number of unanswered questions. For example, what is the maximal size of the subsystem whose Rényi entropy can equilibrate to those of the reduced DE and the reduced CE for a certain initial temperature? And what is the biggest subsystem whose reduced time-evolved density matrix is thermal for a given initial temperature? Why do the eigenvalues of reduced time-evolved states in high-particle-number sub-sectors equilibrate to different values than the reduced DE? Answering these questions needs further detailed study of the diagonal and off-diagonal matrix elements of initial states in the basis of final Hamiltonian and the reduced density matrix of projector operators $|m\rangle\langle m|$ and “off-diagonal” operators $|m\rangle\langle n|$ ($m\neq n$), which can be the subject of future work.

In terms of experimental realization of the models studied, the bosonic model is a simple Bose Hubbard model that can be created with cold atoms on optical lattices. Creating the wall that initially confines the particles to the left half of the system requires local manipulation of the optical potential and may be realized using holographic techniques Islam et al. (2015). The thermalization dynamics from quantum quench for the simple one dimensional Bose Hubbard model with $U/J=0.38$ and effective temperature $\beta J=0.09$ has been studied experimentally Kaufman et al. (2016). The spinless fermions model we studied contains nearest-neighbor and next-nearest-neighbor hoppings and interactions. As discussed in Sec. II, it is equivalent to a zig-zag two-leg ladder with only nearest-neighbor terms (inter and intra-leg). Nearest-neighbor interactions may be created and tuned via Rydberg-dressed potentialsHenkel et al. (2010); Pupillo et al. (2010); Zeiher et al. (2016). The inter and intra-leg nearest-neighbor interactions and hopping can be tuned to different values by varying the lattice spacing on the rung and along the legs to different values and exploiting the range and shape of the Rydberg-dressed potential, as was done in Ref. Zhang et al. (2018). Multi-mode cavity photon-mediated interactions Vaidya et al. (2018) may also be used to create the needed interactions.

Acknowledgements.

This work was supported in part by the National Science Foundation (NSF) under Grant No. DMR-1411345 (SWT) and by the U.S. Department of Energy (DOE) under Award Number DE-SC0019139 (YM).