Thermalization and dynamics in the single-impurity Anderson model

We analyze the process of thermalization, dynamics, and the eigenstate thermalization hypothesis (ETH) for the single-impurity Anderson model, focusing on the Kondo regime. For this we construct the complete eigenbasis of the Hamiltonian using the numerical renormalization group (NRG) method in the language of the matrix product states. It is a peculiarity of the NRG that while the Wilson chain is supposed to describe a macroscopic bath, very few single-particle excitations already suffice to essentially thermalize the impurity system at finite temperature, which amounts to having added a macroscopic amount of energy. Thus, given an initial state of the system such as the ground state together with microscopic excitations, we calculate the spectral function of the quantum impurity using the microcanonical and diagonal ensembles. These spectral functions are compared to the time-averaged spectral function obtained by time evolving the initial state according to the full Hamiltonian, and to the spectral function calculated using the thermal density matrix. By adding or removing particles at a certain Wilson energy shell on top of the ground state, we find qualitative agreement between the resulting spectral functions calculated for different ensembles. This indicates that the system thermalizes in the long-time limit, and can be described by an appropriate statistical-mechanical ensemble. Moreover, by calculating static quantities such as the impurity spectral density at the Fermi level as well as the dot occupancy for energy eigenstates relevant for microcanonical ensemble, we find good support for the ETH. The ultimate mechanism responsible for this effective thermalization within the NRG can be identified as Anderson orthogonality: the more charge that needs to flow to or from infinity after applying a local excitation within the Wilson chain, the more the system looks thermal afterwards at an increased temperature. For the same reason, however, thermalization fails if charge rearrangement after the excitation remains mostly local. In these cases, the different statistical ensembles lead to different results. Their behavior needs to be understood as a microscopic quantum quench only.

Figure 1

Matrix product state (MPS) representation of the kept ($K$) and discarded ($D$) states on the Wilson chain. The $n\mathrm{th}$ box represents the tensor block $X_n$ ($X\in \{K,D\}$), while its bottom, left, and right legs carry the labels of the local states $|{\sigma }_n\rangle$, the kept states ${|se\rangle }_{n-1}^K$, and the kept (discarded) states ${|se\rangle }_n^K$ (${|se\rangle }_n^D$), respectively.

Figure 2

The normalized spectral function $\pi \mathrm{\Gamma }{A}_{\mathcal{E}}$ of the dot level at $T=0$ calculated by using the diagonal ensemble $A_{\mathrm{diag}}$, the microcanonical ensemble $A_{\mathrm{micro}}$, the grand-canonical ensemble $A_{\mathrm{grand}}$, the time-evolved state in the long-time limit $A_{\mathrm{time}}$, as well as the spectral function $A_{\mathrm{\Psi }}$ calculated for the initial state $|\mathrm{\Psi }\rangle =|\mathcal{G}\rangle$. The left inset displays the weights on the Wilson chain $w_n$ for the diagonal ($w_n^{\mathrm{diag}}$), microcanonical ($w_n^{\mathrm{micro}}$), and grand-canonical $\left(w_n^{\mathrm{grand}}\right)$ ensembles for fixed finite length $N=50$, while the right inset represents a zoom of the spectral data around the Kondo peak.

Figure 3

The density matrix weight distributions for the diagonal, microcanonical, and grand-canonical ensembles as a function of the Wilson shell index $n$ [(a), (c), (e)], and the normalized spectral functions [(b), (d), (f)] calculated for a few excited states obtained after acting with the spin-up creation operator of shell $n, {\widehat{f}}_{n\uparrow }^{†}$, on the ground state: (a), (b) $|\mathrm{\Psi }\rangle ={\widehat{f}}_{18\uparrow }^{†}|\mathcal{G}\rangle$, (c), (d) $|\mathrm{\Psi }\rangle ={\widehat{f}}_{23\uparrow }^{†}|\mathcal{G}\rangle$, and (e), (f) $|\mathrm{\Psi }\rangle ={\widehat{f}}_{28\uparrow }^{†}|\mathcal{G}\rangle$. The number of states contributing to the microcanonical ensemble is denoted by $N_{\mathrm{\Psi }}$, where the notation $n=(n_1...n_2)$ indicates the contributing shells. The insets present the zoom of the spectral functions around the Kondo peak.

Figure 4

The matrix elements of the diagonal density matrix ${\sum }_e{\left|C_{nse}^D\right|}^2$, calculated for different excited states $|\mathrm{\Psi }\rangle ={\widehat{f}}_{n\uparrow }^{†}|\mathcal{G}\rangle$, as indicated in the figure panels. The matrix elements are plotted as a function of energy $\omega \equiv E_{ns}^D$ measured relative to the ground-state energy $E_{\mathcal{G}}$. The solid lines display the integrated weights $w_n^{\mathrm{diag}}$ of the diagonal ensemble plotted as a function of energy of a given iteration $n, \omega =\alpha {\mathrm{\Lambda }}^{-n/2}$, with $\alpha$ a numerical constant taking into account the energy spread of each iteration. The vertical dashed lines indicate the average energy of state $|\mathrm{\Psi }\rangle , E_{\mathrm{\Psi }}$.

Figure 5

The normalized spectral function at energy $\omega =0, \pi \mathrm{\Gamma }A\left(0\right)$, calculated using the diagonal and microcanonical ensembles for the initial state $|\mathrm{\Psi }\rangle ={\widehat{f}}_{k\uparrow }^{†}|\mathcal{G}\rangle$ as a function of the Wilson shell index $k$. The spectral data were broadened with ${\omega }_0=E_{\mathrm{\Psi }}$. The predictions based on microcanonical ensemble agree reasonably well with those based on the diagonal ensemble, which indicates that the system behaves thermally.

Figure 6

Energy-resolved expectation values relevant for the microcanonical ensemble of (a) the dot occupation $n_d=\langle E|{\widehat{n}}_d|E\rangle$, and (b) the normalized spectral function operator $\pi \mathrm{\Gamma }\langle E\left|\widehat{\mathcal{A}}\left(0\right)\right|E\rangle$ taken at $\omega =0$. Here, the microcanonical ensemble was based on the initial state $|\mathrm{\Psi }\rangle ={\widehat{f}}_{23\uparrow }^{†}|\mathcal{G}\rangle$, and the energy expectation values are computed for individual energy eigenstates $|E\rangle \equiv {|se\rangle }_n^D$. For (b), the spectral data were broadened using ${\omega }_0=E_{\mathrm{\Psi }}$. The dashed horizontal lines display $\mathrm{Tr}\left[{\widehat{\rho }}_{\mathrm{micro}}{\widehat{n}}_d\right]$ and $A_{\mathrm{micro}}$ for panels (a) and (b), respectively. The left inset in (a) presents the distribution of weights of the microcanonical density matrix $w_n^{\mathrm{micro}}$ on the Wilson chain, where the different symbols differentiate between specific energy shells. The same symbols are also used when plotting the energy expectation values in the main panels. The right insets in (a) and (b) show the corresponding energy-shell-resolved histograms of the data from the main panel.

Figure 7

The weights of the diagonal ensemble $w_n^{\mathrm{diag}}$ as a function of the Wilson shell index $n$ calculated for states $|\mathrm{\Psi }\rangle ={\widehat{n}}_{28\uparrow }|\mathcal{G}\rangle$ and $|\tilde{\mathrm{\Psi }}\rangle =(1-|\mathcal{G}\rangle \langle \mathcal{G}\left|\right)|\mathrm{\Psi }\rangle$.

Figure 8

Weight distribution of various ensembles (upper panels) and the corresponding normalized spectral function $\pi \mathrm{\Gamma }{A}_{\mathcal{E}}$ (lower panels) for the state $|\mathrm{\Psi }\rangle ={\widehat{d}}_{\uparrow }^{†}|\mathcal{G}\rangle$ (left panels) and $|\mathrm{\Psi }\rangle ={|\mathcal{G}\rangle }_{\mathrm{\Gamma }=0}$ (right panels). The grand-canonical spectral function was calculated at an effective temperature $T_{\mathrm{\Psi }}$ specified in the figure, together with ${\omega }_0^{\mathrm{grand}}=T_{\mathrm{\Psi }}$. For the other spectral functions $A_{\mathrm{diag}}, A_{\mathrm{micro}}, A_{\mathrm{time}}$, and $A_{\mathrm{\Psi }}$, we optimized the broadening parameters such that all unphysical features due to discretization were smeared out, with the respective values for ${\omega }_0^{\mathrm{diag}}, {\omega }_0^{\mathrm{micro}}, {\omega }_0^{\mathrm{time}}$, and ${\omega }_0^{\mathrm{\Psi }}$ specified in the panels. The inset to panel (c) shows the weights of the diagonal ensemble $w_n^{\mathrm{diag}}$ plotted on logarithmic scale as a function of Wilson shell index.

Figure 9

(a) The weights of the diagonal density matrix $w_n^{\mathrm{diag}}$ and (b) the contribution to the total energy coming from a given Wilson shell $n, E_{\mathrm{\Psi }}^n=w_n^{\mathrm{diag}}{\omega }_n$, calculated for the state $|\mathrm{\Psi }\rangle ={|\mathcal{G}\rangle }_{\mathrm{\Gamma }=0}$ and plotted as a function of energy $\omega \propto {\alpha }_{\mathrm{\Lambda }}{\mathrm{\Lambda }}^{-n/2}$, with ${\alpha }_{\mathrm{\Lambda }}$ being a constant of the order of unity chosen in such a way that the maximum of weight distribution occurs at the same energy for all values of $\mathrm{\Lambda }$ presented in the figure. The insets show the corresponding data plotted on the $y$-axis linear scale. The vertical dashed lines in (a) and its inset mark the energy scales corresponding to $T_K$ and $\mathrm{\Gamma }$, respectively.

Figure 10

The normalized spectral function for the time-evolved state $A_{\mathrm{time}}$ starting from $|\mathrm{\Psi }\rangle ={|\mathcal{G}\rangle }_{\mathrm{\Gamma }=0}$ calculated at different final times $t_{\mathrm{fin}}$ as indicated in the figure. The spectral data were broadened with a time-dependent broadening parameter ${\omega }_0=\max (1/t_{\mathrm{fin}},T_K)$, with ${\omega }_0=\mathrm{\Gamma }$ for $t_{\mathrm{fin}}=0$, and averaged over time interval $(t_{\mathrm{fin}}-\delta t,t_{\mathrm{fin}})$, with $\delta t=0.2t_{\mathrm{fin}}$. The spectral functions are shifted by 0.01 (0.1) along the $x$ axis ($y$ axis) to increase visibility.

Figure 11

The weight distributions of (a) the diagonal $w_n^{\mathrm{diag}}$, (b) the grand-canonical $w_n^{\mathrm{grand}}$, and (c) the microcanonical $w_n^{\mathrm{micro}}$ density matrices as a function of energy $\omega ={\alpha }_{\mathrm{\Lambda }}{\mathrm{\Lambda }}^{-n/2}$, with ${\alpha }_{\mathrm{\Lambda }}$ of order 1 (see text). These distributions were calculated after application of a single one-particle excitation $|\mathrm{\Psi }\rangle ={\widehat{f}}_{k\uparrow }^{†}|\mathcal{G}\rangle$ at shell $k$, where for a given Wilson shell $k$, the discretization parameter $\mathrm{\Lambda }$ was chosen such that the excitation occurred at comparable energy, i.e., $\mathrm{\Lambda }\left(k\right)=E^{-2/k}$ for $E=T_K$. The vertical dashed lines in (a) present the corresponding energies at which the excitation occurred.

Figure 12

The occupations on the Wilson chain for (a) $\mathrm{\Lambda }=1.5$, (b) $\mathrm{\Lambda }=2$, and (c) $\mathrm{\Lambda }=3$ calculated for initial state (a) $|\mathrm{\Psi }\rangle ={\widehat{f}}_{39\uparrow }^{†}|\mathcal{G}\rangle$, (b) $|\mathrm{\Psi }\rangle ={\widehat{f}}_{23\uparrow }^{†}|\mathcal{G}\rangle$, and (c) $|\mathrm{\Psi }\rangle ={\widehat{f}}_{15\uparrow }^{†}|\mathcal{G}\rangle$ and for diagonal ensemble as a function of the Wilson shell index $n$. The excitation was created at the same energy scale for all $\mathrm{\Lambda }$. The spin-up occupations obtained using ${\widehat{\rho }}_{\mathrm{diag}}$ after having applied $n_x=10$ single-particle excitations are shown for comparison. The inset in (b) presents the buildup of spin-up charge vs shell index $n$ when applying up to a total of $n_x$ excitations. The inset in (c) shows the change of the total charge of the system $\mathrm{\Delta }{N}_{\mathrm{tot}}$ as a function of $n_x$.

Figure 13

The weights of the diagonal $w_n^{\mathrm{diag}}$ density matrix plotted on logarithmic scale as a function of energy $\omega \propto {\mathrm{\Lambda }}^{-n/2}$ calculated for single-particle excitations applied to the ground state ($n_x=1$) and then consecutively to the diagonal density matrix $(n_x>1)$. In (a), ${\widehat{f}}_{k\uparrow }^{†}$ was applied at shell $n_{\max }^{\mathrm{grand}}$, at which $w_n^{\mathrm{grand}}$ is maximum. In (b), the operators ${\widehat{f}}_{k\uparrow }^{†}$ and ${\widehat{f}}_{k\uparrow }$ were applied at shell $n_{\max }^{\mathrm{grand}}$ in an alternating way, whereas in (c) the same operators were applied with probability $p\left(k\right)$ determined by $w_n^{\mathrm{grand}}$. The dashed lines show the corresponding asymptotic lines. After $n_x\gtrsim 5$ excitations the slope of $w_n^{\mathrm{diag}}$ becomes equal to that of $w_n^{\mathrm{grand}}$ (black dashed line). This figure was calculated for $\mathrm{\Lambda }=2.21$, for which the dependence of $w_n^{\mathrm{grand}}$ on energy is $w_n^{\mathrm{grand}}\propto {\omega }^{3.5}$.

Figure 14

(a) Total energy $E_{\mathrm{diag}}$ of the system calculated using diagonal ensemble vs the number of excitations $n_x$ and (b) the weights of the diagonal ensemble $w_n^{\mathrm{diag}}$ as a function of energy after having applied $n_x$ single-particle operators, either ${\widehat{f}}_{k\uparrow }^{†}$ or ${\widehat{f}}_{k\uparrow }$, selected randomly, and applied to shell $k$ with probability $p\left(k\right)$ determined by $w_n^{\mathrm{grand}}$. $T_{\mathrm{grand}}$ is the temperature used for grand-canonical ensemble and $E_{\mathrm{grand}}$ is the corresponding energy of the system. In (a), $E_{\mathrm{diag}}$ for three exemplary runs is shown together with its ensemble-averaged (over 100 runs) value, which shows a steady logarithmic increase with $n_x$ [inset in (a)]. Panel (b) presents the ensemble-averaged weights ${\langle w^{\mathrm{diag}}\left(\omega \right)\rangle }_{100}$ (over 100 runs) of diagonal density matrix plotted as a function of energy $\omega$ with $w^{\mathrm{diag}}\left(\omega \right)\equiv w_n^{\mathrm{diag}}{\mathrm{\Lambda }}^{n/2}$. The dashed line shows the weights of the grand-canonical density matrix $w^{\mathrm{grand}}\left(\omega \right)\equiv w_n^{\mathrm{grand}}{\mathrm{\Lambda }}^{n/2}$ as a function of $\omega$.

Figure 15

MPS representation of the coefficients ${\sum }_e{\left|C_{nse}^D\right|}^2={\sum }_e{\langle \mathrm{\Psi }|se\rangle }_n^D{}_n^D\langle se|\mathrm{\Psi }\rangle$, where the discarded state ${|se\rangle }_n^D$ is defined graphically in Fig. 1. The kept and discarded blocks are denoted by $X=K$ ($X=D$), respectively, while the corresponding blocks of the state $|\mathrm{\Psi }\rangle$ are denoted by ${\psi }_n$ for shell $n$. The star denotes complex conjugation.

Figure 16

The dependence of relative energy uncertainty $\mathrm{\Delta }{E}_{\mathrm{\Psi }}/E_{\mathrm{\Psi }}$ on Wilson shell index for excited states $|\mathrm{\Psi }\rangle ={\widehat{f}}_{k\uparrow }^{†}|\mathcal{G}\rangle$ and for different discretization parameters: (a) $\mathrm{\Lambda }=1.5$, (b) $\mathrm{\Lambda }=2$, (c) $\mathrm{\Lambda }=3$, and (d) $\mathrm{\Lambda }=4$. The horizontal dashed lines indicate the estimate for the relative energy uncertainty $\delta E_{\mathrm{\Psi }}/E_{\mathrm{\Psi }}=1-{\mathrm{\Lambda }}^{-1}$.

Figure 17

MPS representation of a contribution to the time-averaged spectral function of the system due to ${\widehat{d}}_{\sigma }\left(\tau \right){\widehat{d}}_{\sigma }^{†}\left(0\right)$. The open indices need to be summed and multiplied with appropriate exponentials and matrix elements [see Eq. (B2)]. Summation over shells $n=n_0,\cdots ,N$ and states $X{X}^{\prime }{X}^{\prime \prime }\ne KKK$ also needs to be applied.

Figure 18

MPS representation of contributions to element of the new diagonal density matrix ${\rho }_{n^{\prime },s^{\prime }{s}^{\prime }}^{\mathrm{new}}$ [see Eqs. (C4) and (C5)] after having applied a local operator ${\widehat{f}}_{k\sigma }^{†}$ at shell $k$. Note that consistent tracking of environmental degeneracy leads to the additional factor $d^{N-n^{\prime }}/d^{N-n}=1/d^{n-n^{\prime }}$.