Analysis of self-thermalization dynamics in the Bose-Hubbard model by using the pseudoclassical approach

TL;DR

The paper investigates self-thermalization in the Bose-Hubbard model using a pseudoclassical approach to compute the SPDM, revealing that weak interactions induce chaos without destroying BE-like steady-state occupations $n_k=\\frac{1}{e^{\\beta(E_k-\\mu)}-1}$. It shows that two chaotic BH rings equilibrate to a common thermal state and that a chain connected to chaotic reservoirs exhibits a quasi-stationary current consistent with a boundary-driven master equation, with self-thermalization rate $\\\gamma$ scaling as $\\epsilon^2$ relative to a reservoir coupling, placing a bound on Markovian dynamics. The results link microscopic chaotic dynamics to macroscopic thermalization and transport, clarifying time scales for equilibration and the role of reservoir chaos in sustained current. Together, these findings support the applicability of standard statistical mechanics concepts to closed, chaotic many-body BH systems and provide a framework for analyzing thermalization and non-equilibrium transport from first principles.

Abstract

We analyze the self-thermalization dynamics of the $M$-site Bose-Hubbard model in terms of the single-particle density matrix that is calculated by using the pseudoclassical approach. It is shown that a weak inter-particle interaction, which suffices to convert the integrable system of non-interacting bosons into a chaotic system, has a negligible effect on the thermal density matrix given by the Bose-Einstein distribution. This opens the door for equilibration where the two coupled Bose-Hubbard systems, which are initially in different thermal states, relax to the same thermal state. When we couple these two subsystems by using a lattice of the length $L\ll M$, we numerically calculate the quasi-stationary current of Bose particles across the lattice and show that its magnitude is consistent with the solution of the master equation for the boundary driven $L$-site Bose-Hubbard model.

Figures (6)

• Figure 1: Main panel: The line of constant density $\bar{n}=1$, blue dashed curve, and the line of constant energy $\bar{E}=-0.533$, solid read curve. Inset: Bose-Einstein distributions as the function of the quasimomentum $\kappa$ for the points $A$, $C$, and $B$, from bottom to top at $\kappa=0$. Asterisks and open circles show the result of equilibration dynamics for two coupled rings of the size $M=20$, see Sec. \ref{['sec3']}.
• Figure 2: Thermal density matrix for $\beta=2$, $\bar{n}=1$, and $g=0.4$: (a) distribution of 2048 trajectories from the quantum ensemble over the energy; (b) the stationary SPDM in the Bloch basis; (c) Lyapunov exponents of 2048 trajectories; (d) diagonal matrix elements of the stationary SPDM, open circles, as compared to the Bose-Einstein distribution, solid line.
• Figure 3: Equilibration dynamics between two BH systems of the size $M=20$ which are initially in the thermal states with $(\beta,\mu)=(2,-1.07)$ and $(\beta,\mu)=(0.2,-3.63)$, see points $A$ and $B$ in Fig. \ref{['fig00']}. The coupling constant $\epsilon=0.25$.
• Figure 4: Equilibration dynamics for $\epsilon=0.1$, $M=40$, and $L=3$. The interaction constant for bosons in rings is $g=0.4$ while $g=0$ for bosons in the lattice. Initially rings are in the high-temperature thermal states with $\beta_{\rm L}=\beta_{\rm R}=0.2$, and $\bar{n}_{\rm L}=1$ and $\bar{n}_{\rm R}=0.5$. Different panels shows: (a) number of particles in rings and the lattice as the function of time; (b) populations of the individual sites at $t=200$; (c) dynamics of the quantity (\ref{['z']}), red solid line. Additional dashed, dotted, and dash-dotted lines show $z(t)$ for $M=20$ and $\epsilon=0.1/\sqrt{2}$, $\epsilon=0.1$, and $\epsilon=0.1\sqrt{2}$; (d) diagonal elements of the ring SPDMs at $t=0$, solid lines, and $t=200$, markers.
• Figure 5: (a) Lyapunov exponents for $g=0.4$, $\bar{n}=1$, and $M=20$ as the function of the mean kinetic energy and (b) as the function of the total energy; (c) the mean Lyapunov exponent for $g=0.4$, $\bar{n}=1$, and $M=10,20,40,80$, from bottom to top, as the function of the mean kinetic energy (average over 100 initial conditions); (d) Lyapunov exponents for $\bar{n}=1$, $M=20$, and $g=0.2,0.4,0.8$, from bottom to top, as the function of the total energy.